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Abstract 

In  this  thesis,  we  first  study  the  classical  Schwarz  alternating  method  and  an 
additive,  more  parallel  variant  of  it,  known  as  the  additive  Schwarz  method,  ap- 
plied to  solve  saddle  point  linear  systems  obtained  by  discretising  a  saddle  point 
formulation  of  elliptic  Neumann  problems.  We  assume  that  the  discretisation  is 
obtained  by  using  a  mixed  finite  element  method,  in  particular  the  Raviart- Thomas 
elements.  We  prove  convergence  with  a  rate  independent  of  the  mesh  parajneter  h. 
We  aJso  present  results  of  numerical  experiments  using  these  algorithms. 

Following  that,  we  study  two  algorithms  to  solve  problems  on  iteratively  re- 
fined meshes,  namely  the  fast  adaptive  composite  grid  method  (FAC),  and  the 
asynchronous  fast  adaptive  composite  grid  method  (AFAC).  We  give  a  proof  of 
convergence  of  both  these  methods  in  the  mixed  finite  element  case,  with  a  rate  of 
convergence  independent  of  the  mesh  parameters  /i,-.  For  the  FAC  algorithm,  we 
cdso  present  quantitative  bounds  for  the  rate  of  convergence  in  the  case  of  special 
discretisations,  which  show  a  fast  rate  of  convergence,  and  one  which  is  also  inde- 
pendent of  the  geometry  of  the  refined  meshes,  providing  they  axe  shape  regular. 

Finally,  we  study  a  Dirichlet-Neiunann  type  algorithm  for  the  mixed  finite  ele- 
ment case  involving  two  non-overlapping  subdomains.  We  use  this  as  a  precondi- 
tioner  for  a  reduced  Schur  complement  system  obtciined  by  using  an  algorithm  of 
Glowinski  and  Wheeler.  We  prove  that  the  rate  of  convergence  is  independent  of 
h,  the  mesh  parameter.  We  also  present  quantitative  bounds  in  special  geometries, 
that  show  a  fast  rate  of  convergence. 
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Chapter  1 

A  survey  of  background  results, 
old  and  new. 


Introduction. 

In  this  thesis,  we  study  certain  iterative  methods  to  solve  the  saddle  point  linear 
systems  resulting  from  mixed  finite  element  discretisations  of  second  order  linear  el- 
liptic Neumann  problems.  Unlike  standard  finite  element  discretisations  of  elliptic 
problems,  which  lead  to  large,  symmetric,  positive  definite  linear  systems,  the  dis- 
cretisation of  saddle  point  formulations  of  elliptic  problems,  lead  to  large  symmetric 
indefinite  linear  systems.  The  iterative  methods  we  study  are  known  as  domain  de- 
composition methods,  and  are  based  upon  a  division  of  the  domain  of  the  elliptic 
problem  into  smaller  pieces,  so  called  substructures,  with  or  without  overlap.  The 
domain  decomposition  algorithms  involve  the  solution  of  problems  on  each  of  the 
substructures,  often  concurrently,  during  each  iteration.  These  methods  are  therefore 
well  suited  for  parallel  computers.  We  introduce  several  methods,  and  study  their 
convergence  in  the  mixed  finite  element  case. 

In  this  Chapter,  we  present  the  relevent  background  about  saddle  point  formula- 
tions of  elliptic  problems,  and  mixed  finite  element  methods  using  the  Raviart-Thomas 
spaces.  We  then  discuss  an  extension  theorem  for  the  Raviart-Thomas  finite  element 
spaces,  which  we  use  later  in  establishing  bounds  for  the  rate  of  convergence  of  the 
domain  decomposition  methods.  A  section  containing  some  background  on  iterative 
methods,  is  also  included. 

In  Chapter  2,  we  discuss  algorithms  involving  subdomains  with  overlap,  such  as 
the  classical  Schwarz  alternating  method,  cf.  Lions  [22],  and  the  additive  Schwarz 
method,  as  studied  by  Dryja  and  Widlund  [13].    We  present  proofs  of  convergence 


of  these  iterative  methods  when  applied  to  the  mixed  finite  element  case,  and  also 
show  that  the  rate  of  convergence  is  independent  of  the  mesh  parmeter  h.  We  also 
present  numerical  results  of  tests  using  these  methods  with  many  subdomains  and 
a  certain  coarse  mesh  model  problem,  which  improves  the  rate  of  convergence.  The 
results  indicate  a  rate  of  convergence,  which  is  independent  of  the  mesh  parameter  h 
and  even  the  number  of  subdomains.  We  have  also  tested  the  methods  on  problems 
in  which  the  discontinuity  in  the  coefficients  of  the  elliptic  operator  is  large.  Such 
large  variations  in  the  coefficients  occur  in  certain  applications  involving  flow  in  porous 
media.  The  rate  of  convergence  remains  independent  of  the  jump  in  the  discontinmty, 
but  the  accuracy  of  the  pressure  deteriorates.  See  the  section  on  numerical  results, 
in  Chapter  2. 

Following  that,  in  Chapter  3,  we  discuss  certain  iterative  methods,  known  as 
the  fast  adaptive  composite  grid  method  (FAC)  and  the  asynchronous  fast  adaptive 
composite  grid  method  (AFAC),  to  solve  problems  on  repeatedly  refined  meshes. 
These  algorithms  were  originally  developed  for  standard  finite  element  discretisations 
of  elliptic  problems,  cf.  McCormick  and  Thomas  [26],  Mandel  and  McCormick  [24] 
and  Dryja  and  Widlund  [14],  [36].  We  study  these  algorithms  for  mixed  finite  element 
methods.  We  show  that  the  rate  of  convergence  of  the  FAC  algorithm  is  independent 
of  the  mesh  parameters  /i,,  and  also  present  some  quantitative  results  showing  that 
the  rate  of  convergence  is  fast  for  certain  discretisations.  We  also  study  a  variant  of 
the  AFAC  algorithm,  and  show  that  it  has  a  rate  of  convergence  independent  of  the 
mesh  parameters  h,. 

Finally,  in  Chapter  4,  we  study  a  domain  decomposition  method  involving  two 
subdomains  without  overlap,  in  the  mixed  finite  element  case.  This  is  based  on 
an  algorithm  of  Glowinski  and  Wheeler  [16],  and  ideas  about  a  Neumann- Dirichlet 
preconditioner  used  in  standard  finite  element  discretisations  of  elliptic  problems,  cf. 
Bj0rstad  and  Widlund  [4]  and  Bramble,  Pasciak  and  Schatz  [5].  See  also  Quarteroni 
[29]  for  related  work  on  saddle  point  problems.  We  prove  that  the  Dirichlet- Neumann 
algorithm  has  a  rate  of  convergence  independent  of  the  mesh  parameter  /i,  and  also 
indicate  that  for  certain  geometries  a  fast  rate  of  convergence  is  obtained. 


1.1      Function  spaces  and  trace  theorems. 

In  this  section,  we  introduce  some  function  spaces  and  also  present  trace  and  exten- 
sion theorems  which  will  be  used  later  on.  All  of  the  spaces  introduced  are  Hilbert 
spaces.  First  we  introduce  the  Sobolev  spaces  and  some  trace  and  extension  theo- 
rems; cf.  Lions  and  Magenes  [21],  Necas  [27],  and  Grisvard  [19].  Then  we  introduce 
the  H{div,.)  space,  some  of  its  subspaces,  and  related  trace  and  extension  theorems. 
See  Girault  and  Raviart  [15],  and  Raviart  and  Thomas  [30]. 

1.1.1      Sobolev  spaces. 

Let  Cl  C  R^  he  a.  Lipschitz  region,  i.e.,  dO.  can  be  represented  locally  as  a  Lipschitz 
function.  We  denote  by 


H°{n)  =  L\n)  =  {u:  I  \u\''dx  <  oo}, 


and  define  the  norm  by  ||u||L2(n)  =  {JnU^dxY^^.  For  integer  m  >  0,  we  define  the 
Sobolev  space  by, 

H'^iQ)  =  {ii  e  L\n)  :  d°u  G  L\n),  for  |a|  <  m}, 

where  a  =  (qi,Q2)  are  non-negative  integer  indices  representing  the  order  of  the 
partial  derivitives;  |q|  =  ai  +  Q2;  and  d'^u  =  {■£-^)°''{-^)°^u.  The  norm  and  semi- 
norm  are  defined  respectively  by 


1/2 

ia|<m  '  "  |a|  =  m 

For  non  integer  s,  with  m  <  5  <  m  +  1,  where  m  is  a  non-negative  integer,  we  can 
define  the  fractional  order  Sobolev  space  H'{Q,),  and  its  associated  norm  by  using  the 
theory  of  interpolation  spaces, 

These  spaces  form  a  Hilbert  scale,  and  satisfy  various  interpolation  inequalities,  see 
Lions  and  Magenes  [21]. 

For  0  <  e  <  1,  we  also  have  the  following  equivalent  definition  of  the  fractional 
order  Sobolev  space  H'{Cl): 

H'iQ)  =  {ue  L\n)  :  |u|„.(n)  <  00}, 

3 


equipped  with  an  equivalent  semi-norm  and  norm  defined  respectively  by 
kUno)     =    {LenJyeuiH^)-<y)\'/k-yr'')dxdyy'', 

For  0  <  5  =  m  +  e,  where  m  >  0  is  an  integer  and  0  <  e  <  1,  we  define  the  fractional 
order  Sobolev  space  H'{Cl)  by, 

H'{n)  =  {ue  H"'{n)  :  d^'u  e  H'{n),     \a\  =  m}, 
and  define  an  equivalent  norm  by 

||'"llH'(n)  =  (ll'"llH'"(n)+    H   Id^ulH.^n))^^^. 

\a\  =  m 

Having  defined  the  Sobolev  spaces  H'{Q.)  for  non-negative  s,  we  now  present  the 
definition  of  their  dual  spaces.  For  5  >  0,  we  define  the  dual  space  {H'{Q,))'  as  the 
space  of  continuous  linear  functionals  on  H'(Q): 

iH'in))'  =  {v{.)  :  H'{n)  -^  R;\\v{.)\\  <  oc}, 

with  the  norm  defined  by 

II    II  _  IJuUvdxl 

\\v\\{H'{n)y  =     sup 


where  /q  uvdx  represents  the  duality  pairing.  These  spaces  form  a  Hilbert  scale  with 
a  dense,  compact  imbedding: 

for  5i   <  52- 

Next,  we  discuss  the  Sobolev  spaces  of  functions  on  the  boundary  dCl.  Since,  we 

will  need  to  use  these  boundary  spaces  only  for  —  1  <  5  <  1,  we  present  their  definition 

only  for  these  vaJues.  L^{dQ.)  is  defined  using  the  standard  definition  of  the  boundary 

integrcd: 

L^dn)  =  {u:  I     \u\^ds^  <  oo}, 
Jan 

with  the  norm  defined  by 

\\uU^ao)  =  {[    u'ds^Y'\ 
Jan 


For  0  <  5  <  1  the  Sobolev  spaces  on  the  boundary  is  defined  by 

H'{dn)  =  {ix  e  L\dn)  -.  |-u|H'(an)  <  oc},  (1.1) 

where  the  semi-norm  and  norm  are  defined  respectively  by 

WlH'iau)    =    iUaoJyeaoiH^)-'^{y)\V\x-yr'')ds^dsyY/\ 

(1.2) 

For  —1/2  <  5  <  0,  we  define  the  Sobolev  space  H'{dVL)  by  duality. 

Let  r  C  dVt  be  a  nonempty  proper  subset  of  dVt.  For  s  =  |,  we  introduce  the 
interpolation  space 

Hli\V)^[L\V),nl{V)]y,. 

The  space  HIq\V)  f  H^'\l),  but  is  continuously  imbedded  in  H^I^{V): 

and  has  the  property  that  functions  in  Hqq  (F)  can  be  extended  by  zero  to  the  rest 
of  dQ,  as  a  continuous  map,  to  a  function  in  n^^^{dQ).  An  equivalent  definition  of 
the  .ffoo  (T)  norm  is  obtained  by: 

where  t  is  an  arclength  parameter  along  F  and  L  is  the  length  of  F.  One  can  verify 
using  this  formula  and  the  definition  of  the  H^^^{dQ.)  norm  given  by  equation  (1.1), 
that  a  function  in  H^^^{dQ)  which  is  zero  on  dO.  -  F,  has  its  H^^^{dn)  norm  bounded 
by  the  HII)^{T)  norm.  Using  this,  we  also  see  that  functions  in  H'^^^{T),  when  extended 
by  zero  to  the  rest  of  dQ,  do  not  neccessarily  belong  to  H^^^{dQ.).  The  dual  space  will 
be  denoted  by  {Hll^{T)y.  See  Lions  and  Magenes  [21],  for  details  about  the  fl'oo^(.) 
spaces. 

Next,  we  define  the  trace  map,  and  state  a  trace  theorem.  For  smooth  functions 
defined  on  fi,  we  have  an  obvious  pointwise  definition  of  its  trace  or  boundary  values 
on  5J1.  The  following  Lemma  extends  this  notion  of  trace  to  H'{Cl)  functions.  Since 
n  C  il^  is  a  Lipschitz  region,  we  state  the  result  only  for  1/2  <  s  <  3/2,  however  we 
remark  that  for  smooth  domains,  the  result  is  true  for  all  5  >  1/2.  For  a  proof,  see 
Necas  [27],  or  Grisvard  [19]. 


Lemma  1  (Trace)   Let  1/2  <  s  <  3/2.     Then  the  trace  map  7,  with  ^u  =  u  |en 
for  smooth  functions  u,    can   be    extended  as   a   continuous   map  from  H'{Cl)    onto 

H'-'l\d^) 

with  a  positive  constant  c{Q,s)  such  that 

||7u||;^.-i/2(an)  <  c{n,s)\\u\\H'{Q)- 
The  constant  c{Q,s)  depends  upon  s,  and  the  Lipschitz  constant  of  the  region  J). 

Since  7  maps  H'{Q.)  continuously  onto  F'-^/^(5fi),  we  can  apply  the  open  mapping 
theorem  to  get  a  continuous  right  inverse  E 

E:H'-'l\d^) — >  H'[^),  (1.3) 

with 

lEg  =  g,    yg  e  H'-'/'{dn). 

This  is  an  extension  theorem. 

Using  the  trace  theorem,  we  are  able  to  define  a  closed  subspace  of  H^{D,).   We 
let  H^{n)  denote 

H'o{Q)  =  {ue  H\n):-fu  =  0  on  on}, 

and  it  is  provided  with  the  H^{^)  norm.  We  now  state  Freidrichs'  inequcJity  about 
functions  in  Hq{CI): 

Lemma  2  (Freidrichs*  inequality)    There  exists  a  positive  constant  c(Q)  such  that 
for  all  ue  H^{n): 

|l^||H>(n)  ^  c(Q)|t/|Hi(n)- 

We  also  present  Poincare's  inequality. 

Lemma  3  (Poincare's  inequality)    Tlicre  exists  a  positive  constant  c(Cl,s)  such 
thai  for  all  u  6  H'{VL),  for  1  >  s  >  0,  we  have 


\\u\\H'W<c{n,s)[\u\],.^^^  +  {f  udxYY^' 


Let  Vr{^)  denote  the  polynomials  of  degree  less  than  or  equal  to  r  on  fi.  Then,  by 
us'ng  Poincare's  inequality,  it  can  be  shown  that  the  fi^'-seminorm  is  equivaler ".  to 
the  ^'-norm  on  the  quotient  space  H'{Q,)l'Po{Vt).  For  a  proof  of  this  result,  see  Necas 
[27]. 

We  end  this  section,  by  introducing  the  Sobolev  spaces  (J?*(fi))'  of  vector  func- 
tions, u.  Let  u  =  {ui,U2).  We  associate  the  following  norm  on  {E'{Q)y  for  5  >  0: 

\m(H'{n)y  =  (||'"i||^.(n)  +  IMH-^n))^^^- 
The  dual  spaces  are  defined  analogously  to  the  scalar  case. 

1.1.2      The  H{div,Q)  space  and  its  subspaces. 

We  now  introduce  a  function  space  that  will  be  used  in  the  mixed  formulation  of 
elliptic  problems.  Let  H{div,Q)  denote 

equipped  with  the  following  norm 

Clearly,  {H^{Q)Y  is  continuously  imbedded  in  H{div,Q). 

We  now  discuss  the  normal  component  of  a  vector  function  on  the  boundary  of  its 
domain  of  definition.  For  a  Lipschitz  domain  Q,  the  unit  normal  n  to  the  boundary 
d^  exists  almost  everywhere.  Thus,  for  a  smooth  vector  function  u  on  Q,  the  normal 
component  u-n  on  dO.  exists  almost  everywhere.  The  following  normal  trace  Lemma 
extends  the  notion  of  the  normal  component  to  H{div,Cl)  functions.  We  include  a 
proof.  Also,  see  Raviart  and  Thomas  [30]. 

Lemma  4    Given  u  G  H{div,Cl),  there  exists  a  map  fn, 

fr.:H{div,n)^H-'^'{dn), 

such  that  for  smooth  vector  functions  u,  we  have  7„u  =  u  ■  n  on  dO,.  There  exists  a 
constant  c(Q),  which  depends  on  the  Lipschitz  constant  of  the  domain  Q,,  such  that, 

||7r.^ilH->/=(an)  <  c(f))|]u||^(j,„n),      Wu  e  n{div,n). 


PROOF  OF  THEOREM.  Let  g  G  H^^^{dQ,).  Let  (f)  =  Eg,  be  its  extension  into  F^(fi),  as 
defined  by  equation  (1.3).  Then  , 

V  •  {4>u)  =  (p{V  ■u)^V(p-u. 

Integrating  and  using  the  divergence  theorem,  we  obtain 

/     g^^udsj:  =       {4)V  •u  +  V(})-u)dx, 
Jau  Jn 

since  n  ■  u  =  7„it.  Applying  the  Cauchy-Schwarz  inequality,  we  obtain 

I  /    g-inuds^]  <  \\(i>\\m[U)\\u\\ji{div,ny 

JdU  ^  ' 

and  using  the  the  boundedness  of  the  extension  map  E,  we  obtain 

Using  the  definition  of  the  H'^^'^{dQ)  norm,  we  obtain 

||7„u|i^-./2(3n)  <  4^)\MH{div,Q)- 

This  proves  that  7^  is  continuous. 

We  now  show  that  7n  maps  H{div,Q)  onto  H~^^^{dQ).  To  do  this,  we  pose  the 
following  Neumann  problem  for  (^  with  Neumann  data  g  G  H~^^^{dQ),  using  its 
variational  formulation  on  H^{Q,): 

(   _A(5i  +  (^    =     0      in  n, 
\        d4>ldn     =    g      \xi  dO., 

By  construction,  4>  £  H^{Vl).  Also,  from  the  equation  we  have  V(j)  G  H{div,Q,)  with 
7„V*(/.  =  g.  Thus  7„  maps  H{div,n)  onto  H-^/^{dn).  □ 

By  an  application  of  the  open  mapping  theorem,  we  have  a  continuous  right  inverse 
En,  extending  H~^^^{dQ,)  functions  into  H{div,Q). 

Next,  we  define  the  divergence  free  subspaces  of  the  H{div,Q,)  spaces: 

H{div°, n)  =  {ue  H{div,  Q)  :  V  •  u  =  0}. 

H{div°,Cl)  is  equipped  with  the  H{div,Q)  norm.  We  note  that  the  H{div,Cl)  norm 
is  equal  to  the  (i^(f2))^  norm  on  H{div°,Q). 


Next,  we  define  subspaces  of  H{div,U)  with  zero  normal  trace  on  dO,  : 

HofiQidiv,  Vt)  =  {u  E  H{div,  Q)  :  7„u  =  0}, 

equipped  with  the  H{div,Q,)  norm.  Similarly,  we  denote  by  no,ao{div°,Q): 

HoM<ii^°,^)  =  {ue  H{div\n)  :  7„u  =  0}, 

equipped  with  the  H{div,Q)  or  {L'^{Q)Y  norm. 

Finally,  we  state  a  result  about  composite  functions.  See  Raviart  and  Thomas  [30]. 
Let  Q  =  U^if2j  where  the  {f^j}  are  mutually  disjoint.  Then  we  have  the  following 
neccessary  and  sufficient  condition  for  v  to  be  in  H{div,Cl): 

Lemma  5   Suppose  that  f  |n,  G  H{div,  fij)  for  i  =  1, . . . ,  N ;  and  suppose  that 

f  n,j  ■  v\n,  +  Uj,  ■  v\n,  =  0  on  dil,  n  dQj     \/i,j,  ,^  ^. 

1    where  ri,j  is  the  outward  normal  from  Q,  to  fl^  on  5fi,  H  dClj. 

Then  V  e  H{div,n). 

PROOF  OF  LEMMA .  Let  v  satisfy  the  above  hypothesis.  We  will  show  that  v  G  H{div,  Q). 
First,  since 

v\a,e{L'{n,))\yi^ve{L'in))\ 

Next,  we  will  show  that  the  divergence  of  v,  in  the  sense  of  distributions,  is  in  L^{Q,), 
thereby  making  v  G  H{div,Q).  Our  candidate  for  V  •  v,  in  the  sense  of  distributions, 
is,  simply  the  L^{Q.)  function  d{x)  obtained  by  letting: 

d{x)\u,  =  v-v{x)\a„    eL\n,), 

for  each  fi,.  Let  </>  G  C^{0.)  be  an  infinitely  differ enti able  test  function  with  compact 
support  in  Q.  Then,  by  the  definition  of  the  distributional  derivitive,  we  have 

/  (pV  •  vdx  =  —       V<f)  •  vdx. 

Thus,  we  need  to  show  that 

I   (f)V  -vdx-  I   (f)d{x)dx  =  0. 


Using  the  definition  of  d(a;)|n,  and  of  the  distributional  derivitive  and  integrating  by 

parts,  we  obtain, 

N     . 

f  4>V  •  vdx  =  -  /  W  ■  ^dx  -  -  I   V4)-  vdx  -  X!  /  ^^  ■  ^dx, 
Jn  Jn  Jn  ^-j  •'n, 

=  -       ys  .  vdx  +  y^  /  V0  •  vdx  -Y\  I      4>^  ■  vdx 
Jn  ~{ Jo,  f^^  Jaa, 

=  —  y^  /  (hn  •  vdx. 

Now  by  assumption,  we  have  a  cancellation  of  the  sum  of  normal  traces  of  v  on 
dVl,  n  d^j.  Thus,  this  last  term  is  zero.  This  proves  that  V  ■  v  =  d{x),  and  therefore, 

V  • -^  G  X==(n),  since  cf(x)  e  ^^(Q).  Thus,  v  e  ^(Jzr,  f^). 

The  neccessary  part  of  the  theorem  is  easily  proved  by  using  the  normal  trace 
Lemma  ,  and  the  property  that  the  restriction  of  a  H{div,  Q)  function  to  a  subdomain 
Q,  is  in  H{div,Cl,).  □ 

1.2      A   saddle-point   formulation  of  elliptic  prob- 
lems. 

In  this  section,  we  introduce  the  elliptic  problem  that  we  consider  throughout  this 
thesis.  First  we  describe  the  elliptic  partial  differential  equation.  We  also  present  an 
abstract  framework  which  is  used  to  obtain  the  existence  and  uniqueness  of  solutions 
of  the  saddle  point  problem. 

Then  we  describe  its  saddle-point  variational  formulation.  Next,  we  describe  how 
a  discrete  approximation  to  the  saddle-point  problem  can  be  obtained  by  restricting 
the  variational  problem  to  finite  dimensional  subspaces.  Finally,  we  describe  the  par- 
ticular finite  dimensioned  spaces  we  use,  namely,  the  Raviart-Thomas  finite  element 
spaces. 

Consider  the  following  elliptic  problem  for  p  in  a  2  dimensional  polygoncd  region 

I  -V-MVp)    =    /     inn 

\       n-{AVp)     =    g      in5Q,  ^       ' 

where  /  and  g  satisfy  the  compatability  condition  /^  fdx  +  JQ^gds^  =  0,  and  A{x) 

is  a  2  X  2  symmetric  matrix  function  with  i°°(n)  components  satisfying 

Am,n(^(a:))  >  Q  >  0, 
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u    =    AVj> 

in  CI 

V-u    =    f 

in  n 

n-u     =     g 

in  dn 

almost  everywhere,  for  some  positive  constant  a. 

Note  that  specifying  the  conormal  derivitive  of  p,  n  •  AVp  =  g,  on  dCl,  is  equal  to 
specifying  the  normal  derivitive,  dp/dn  =  g  on  dQ,  when  A{x)  =  I,  the  identity.  In 
applications  to  fluid  flow  in  porous  medium,  where  p  denotes  the  pressure  and  AVp 
denotes  the  fluid  velocity  given  by  Darcy's  law,  the  conormal  derivitive  of  p  is  the 
same  as  the  normal  trace  of  the  velocity  on  dCl. 

In  such  applications,  approximations  to  the  fluid  velocity  u  can  be  obtained  by 
discretising  the  elliptic  problem  in  p  and  then  A{x)'Vp.  Alternatively,  we  can  form  an 
equivalent  system  using  both  p  and  AVp  as  unknowns  and  study  discretisations  of  the 
resulting  system  of  partial  differential  equations.  We  present  some  details.  Introduce 
u  =  AVp  as  an  additional  unknown.  Then  we  obtain  the  transformed  system: 


(1.6) 


We  may  use  finite  difference  formulas  to  discretise  the  enlarged  system  (1.6).  However, 
in  this  thesis  we  consider  only  discretisations  obtained  by  using  the  Raviart-Thomas 
finite  elements,  which  in  some  cases  are  equivalent  to  finite  difference  approximations, 
see  Wheeler  and  Gonzalez  [33].  We  introduce  the  variational  problem  associated  with 
system  (1.6).  Without  loss  of  generality,  we  assume  that  ^  =  0.  Multiplying  equation 
(1.6)  by  suitable  test  functions,  and  integrating  by  parts,  we  obtain  the  following 
equivalent  variational  form: 

Find  u  e  Ho,dCt(div,  Q),  p  €  L^{Q.)/R  such  that 

J^u^A-\x)vdx  +  J^p{V  ■v)dx     =     0,  W  e  Ho,an{div,n)         (1.7) 

J^qiV-u)dx  =     -Sofqdx,      VqeL'{Cl)/R. 

We  now  discuss  an  abstract  framework  for  studying  the  saddle-point  problem 

(1.7).  See  Giraiilt  and  Raviart  [15],  for  more  details.  Let  X  and  Y  be  Hilbert  spaces, 

let  a(.,  .),fc(., .)  be  continuous  bilinear  forms,  and  let  W{.),F{.)  be  continuous  linear 

functionals,  defined  on: 

a:  X  xX  — y  R 

b:  X  xY  — >  R 

W  :  X  —^  R 

F  :  Y — y  R. 

We  then  have  the  following  existence  and  uniqueness  theorem  for  the  following  prob- 
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lem: 

Find  u  G  A',  p  e  y  such  that 

a{u,v)  +  b{v,p)  =   w{v),     yvex  (1.8) 

b{u,q)  =     F{q),       VgGy  , 

Theorem  1   Let  a,b,W,F  be  defined  as  above.    Suppose  there  exists  a  positive  con- 
stant a  such  that 

a{u,u)>a\\u\\\,     yueXo,  (1.9) 

i.e.,  a(.,.)  is  Xq- elliptic,  where 

Xo  =  {veX:b{v,q)  =  0,    Vg  G  Y}.  (1.10) 

Also  suppose  that  there  exists  a  positive  constant  /?  such  that 

v^ey,    sup^>/3||g|ly,  (1.11) 

v€A'    \\V\\X 

I.e.,  b{.,.)  satisfies  the  Babuska-Brcz:i  \n{-sup  condition.    Then  problem  (1.8)  is  well 
posed,  I.e.,  has  a  unique  solution  {u,p)  £  A'  X  Y ,  satisfying  the  following  bounds: 

\\u\\x<c{a,(3){\\F\\x>  +  \\W\\Y.), 
\\p\\Y<c{a.f3){\\F\\x-+\\W\\Y-). 

The  constant  c(q,/3)  depends  only  on  a  and  (3. 

REMARK:  By  the  continuity  of  the  bilinear  forms  a(., .)  and  6(., .),  there  exists  contin- 
uous linear  maps 

A:X  V  X' 

B  -x  — >  y, 

defined  by 

{Au.v)  =  a{u.,v).,      Vu,r  G  A, 
{Bu,q)  =  b{u,q),       WueX,qeY, 

where  (., .)  denotes  the  action  of  a  linear  functional. 

REMARK :  The  inf-sup  condition  (1.11)  is  equivalent  to  the  map  B  being  an  isomorphism 

between 

B-.X^—^Y', 

with 

\\Bv\\y'>(3\\v\\x,     VvGXo^ 
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Thus  if  the  inf-sup  condition  holds,  then  B  maps  X  onto  Y'  ajid  has  a  continuous 
right  inverse,  the  norm  of  which  can  be  estimated  by  1/(3.  Thus,  given  q'  €  Y\  there 
exists  V  e  X  such  that  Bv  =  q'  satisfying 

Mx  <  {i/(3W\\y- 

REMARK:  The  inf-sup  condition  (1.11)  is  also  equivalent  to  the  map  B^  being  an 
isomorphism  between 

B^:Y—^{X^y, 

with 

ll5^g|U'>/3|kl|y,     Vg  G  Y. 

Thus,  if  the  inf-sup  condition  holds,  and  u'  G  {X^)',  i.e.,  {u',v)  =  0,  for  all  v  €  Xo, 
there  exists  p  e  Y  with  B^p  —  u'  and  satisfying 

IIpIIv  <  (l//3)|lu'lU.. 

To  apply  this  abstract  framework  to  the  mixed  formulation  of  the  elliptic  problem 
(1.7),  to  obtain  well  posedness,  we  need  to  check  the  Xo-ellipticity  and  the  inf-sup 
condition  for  problem  (1.7).  In  our  application,  we  have 

X  -  Ho,da{div,Q) 

Y  =  L^{n)/R 

a{u,v)  =  J(^v{xfA-\x)u{x)dx  ^_^2) 

b{u,p)  =  f^{V  ■  u{x))p{x)dx 

W{v)  =  0 

F{q)  =  -Jnfix)qix)dx, 

and, 

Xo  =  Ho,an{div°,^). 

First  we  check  that  a(.,.)  is  Xo-elliptic.  By  assumption,  \„,in{A{x))  >  a  >  0,  for 
almost  every  x.  Therefore,  since  for  divergence  free  functions,  the  H{div,.)  norm  is 
equal  to  the  (i^(.))^  norm,  we  have 

a(^,u)>a||u||i(,,^^^^,     yueH{div°,n).  (1.13) 

Thus  a(.,.)is  Ho.da{div°,n)-empt\c. 
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Next,  we  check  the  inf-sup  condition.  A  Lemma  in  Girault  and  Raviart  [15],  states 
that  given  q  G  L^{Cl)/R,  there  exist?  u  G  (i/o(^))^  "^i^h  V  •  u  =  g  and  satisfying 

||^||(H'(n))2  <  c||g||LJ(n), 

for  a  positive  constant  c.   By  the  imbedding  of  {H^{Q.)y  in  H{div,D,),  we  have  u  6 
H{div,Q)  and  we  obtain  an  inf-sup  condition  for  problem  (1.7). 

Thus  the  variational  problem  (1.7),  the  mixed  formulation  of  an  elliptic  Neumann 
problem,  has  a  unique  solution  in  Ho,drt{div,^)  x  L^{Cl)/R.  Note,  since  we  are  study- 
ing a  Neumann  boundary  value  problem,  the  pressure  is  unique  only  in  L^{Cl)/R,  i.e. 
unique  up  to  a  constant  in  L^{Q). 

1.3      An  abstract  framework  for  the  discretisation 
of  the  saddle-point  problem. 

Now  that  we  have  introduced  the  abstract  variational  problem  (1.8),  and  the  specific 

problem  (1.7),  that  we  study  in  this  thesis,  we  discuss  how  to  discretise  the  problem. 

For  each  h,  let  A';,  and  Yh  be  finite  dimensional  subspaces  of  X  and  Y  respectively. 

For  each  h  we  approximate  the  variational  problem  by  restricting  it  to  Xh  x  Yh  as 

follows: 

Find  Uh  G  A'^,  ph  G  Yh  such  that 

<     a{uh,v)-rb{v,ph)  =     W{v),      ^v  e  Xh  (1-14) 

,   b{uh,q)  -     F{q),        yqeYh 

Thus,  for  each  /i,  (uhiPh)  will  be  considered  as  an  approximation  to  {u,p),  the  solution 

of  problem  (1.8).  Here  h  denotes  a  parameter,  usually  the  maximum  mesh  width, 

which  tends  to  zero. 

For  the  discrete  variational  problem  (1.14)  to  be  well  posed,  we  assume  that  for 

each  h  there  exists  positive  constants  ah^^h-,  possibly  dependent  on  /i,  such  that 

the  bilinear  forms  a(., .)  and  6(., .)  satisfy  the  AT/i^o-ellipticity  and  the  Babuska-Brezzi 

inf-sup  condition  respectively: 

\/vheXhfi,      o-{vh-,Vh)  >  ah\\vh\\]c, 

^qeYh,         sup  ^[^  >  ^,||5|iy, 
vex,  \\v\\x 

where 

Xhfi  =  {vhE  Xh:b{vh,q)  =  0,\/q  eYh}.  i 
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Then  we  have  well  posedness  of  problem  (1.14)  by  Theorem  1.  In  addition,  see  Girault 
and  Raviart  [15],  we  have  the  following  error  bound: 

\\u  -  Uh\\x  +  \\p  -  PhlW  <  c{ah,0h){  inf    \\u  -  Vh\\x  +   inf   ||p  -  9/.||y}        (1.15) 

where  {u,p)  are  solutions  of  the  continuous  problem  and  {uh,Ph)  are  solutions  of 
the  discrete  problem.  In  most  applications,  if  the  spaces  Xh  and  V/,  are  chosen  to 
approximate  well  X  and  Y  respectively,  then  the  approximation  error  can  be  shown 
to  approach  zero  as  h  approaches  zero. 

This  discretization  leads  to  a  saddle  point  problem,  i.e.,  a  symmetric  indefinite 
linear  system  with  the  following  block  structure: 

(^:f)(:)=('^:)- 

Ah  is  symmetric  and  positive  definite  and  approximates  the  map  A.  It  is  uniformly 
positive  definite,  (in  h),  in  the  subspace  {vh  ■  BhVh  =  0}.  B^  is  the  discrete  ap- 
proximation to  the  B  map.  Thus  {uh,Ph)  can  be  found  by  solving  the  linear  system 
(1.16). 

We  mention  here  that  the  spaces  A';,  and  Yh  should  be  chosen  so  that  the  compu- 
tation of  the  stiffness  matrices  (1-16)  is  tractable  and  further,  that  they  approximate 
the  original  problem  well.  One  such  choice  is  the  use  of  finite  element  spaces.  Here 
we  describe  briefly  their  construction.  Let  us  assume  that  the  spaces  X  and  Y  in 
the  abstract  variational  framework  of  Theorem  1  are  function  spaces  defined  on  a 
polygonal  domain  J7.  For  each  h,  where  h  is  a.  real  parameter,  we  partition  the  do- 
main into  triangles  (quadrilaterals)  with  the  width  of  each  triangle  (quadrilateral) 
being  less  than  h.  Then  we  define  the  spaces  Xh  and  Yh  as  consisting  of  piecewise 
smooth  functions,  usually  polynomial  in  each  triangle  (quadrilateral),  and  satisfying 
some  global  regularity.  In  particular,  throughout  this  thesis,  for  the  mixed  formu- 
lation of  elliptic  problems,  (1.7),  we  use  the  Ravi  art -Thomas  finite  element  spaces. 
For  r  =  0, 1,2,3, . . .,  these  finite  element  spaces  contain  polynomicJs  of  order  r,  and 
satisfy  a  uniform  inf-sup  condition,  i.e,  these  spaces  satisfy  the  inf-sup  condition  with 
a  positive  constant  j3h  independent  of  h.  We  discuss  these  finite  element  spaces  in 
more  detail  in  the  next  section. 
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1.4      The  Raviart- Thomas  finite  element  spaces. 

In  this  section,  we  introduce  the  Raviart-Thomas  finite  element  spaces,  see  Raviart 
and  Thomas  [30],  which  we  use  to  discretise  the  mixed  formulation  of  elliptic  problems 
(1.7).  Throughout  the  rest  of  the  thesis,  we  assume  that  fi  is  a  polygonal  region  in  R^. 
We  partition  Q  into  disjoint  subregions  called  elements,  which  are  either  all  triangles 
or  all  parallelograms.  Let  /i  be  a  parameter  denoting  the  maximum  diameter  of 
elements  associated  with  the  partition.  We  denote  the  partition  or  triangulation  by 
r'',  and  the  individual  elements  by  K.  The  Raviart-Thomas  finite  element  spaces, 
Vh{Cl)  and  Qh{^),  based  on  the  triangulation  t^,  are  subspaces  of  H{div,Q)  and 
X^(n),  respectively,  consisting  of  functions  which  are  polynomial  in  each  element 
K,  and  satisfy  certain  conditions.  We  describe  them  in  more  detail  in  the  following 
subsections.  We  discuss  two  types  of  spaces,  one  based  on  triangular  elements,  and  the 
other  based  on  quadrilateral  elements.  We  first  define  the  spaces  V/j(0)  and  Qhi^), 
which  are  valid  for  Dirichlet  boundary  conditions.  Then,  we  choose  subspaces  Xh{^) 
and  Yh{^)  of  Vh{Cl)  and  Qhi^)  respectively,  which  are  valid  for  Neumann  boundary 
value  problems,  the  specific  problem  we  consider  throughout  this  thesis. 

1.4.1      Triangular  elements. 

As  before,  we  let  f]  be  a  polygonal  region  in  R^ .  We  partition  it  into  triangles 

where  t^  denotes  the  triangulation  of  Q,  as  before.  We  assume  that  the  vertices  of  the 
triangles  either  coincide  with  the  vertices  of  other  triangles  or  lie  on  the  boundary  dCl; 
In  particular,  we  assume  that  the  vertices  do  not  lie  on  the  edges  of  other  triangles. 
For  each  triangle  K,  we  denote  by 

hfi  =    Diameter  of  K, 

pji  =    Diameter  of  circle  inscribed  in  K . 

DEFINITION .  A  family  of  triangulations  {r''}  of  f2  is  said  to  be  shapt  regular^  if  there 
exists  a  constant  ai  >  0,  such  that, 
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DEFINITION.     A  family  of  triangulations  {t^}  of  CI  is  said  to  be  uniformly  shape 
regular,  if  there  exists  positive  constants  (Ti,«r2,  such  that: 

\/h,\/K^T^,        (TiPK  >  fiK  >  <T2h. 

Affine  maps  and  correspondences. 

Let  K  denote  the  unit  triangle  with  vertices 

ai  =  (0,0);    a2  =  (l,0);    03  =  (0, 1). 

Given  a  non-degenerate  triangle  K  E  t^  with  vertices  01,02,03,  we  map  K  onto  K, 
so  that  FA'(a,)  =  a,,  for  i  =  1,2,3,  using  an  invertible  affine  linear  map  Fk'- 

of  the  form  Fxii)  =  B/ci  +  ^k,  where  Bk  is  a  2  x  2  invertible  matrix  and  bjc  is  a 
2-vector.  Using  the  map  Fk(.),  we  can  define  a  one-to-one  correspondence  between 
scalar  and  vector  functions  on  K  and  K. 

For  any  scalar  function  4>  defined  on  ii',  we  associate  the  scalar  function  4>  defined 
on  K  by  a  one-to-one  correspondence  <p  < >•  (p,  given  by: 

4>  =  4>oFji\     {4>  =  <i>oFK). 

For  vector  functions  v  =  {vi^v^),  defined  on  iT,  however,  to  preserve  the  normal  trace 

of  certain  functions,  we  define  v  on  K  \>y  s.  one-to-one  correspondence  v  < >  v  given 

by: 

V  =  {1/Jk)Bkv  0  Fj^\      {i  =  JkBJ,'vo  Fk), 

where  Jk  =  det{BK)- 

The  following  Lemma,  found  in  Raviart  and  Thomas  [30],  states  some  of  the 
properties  of  the  divergence  and  the  normal  trace  of  vector  functions  under  the  cor- 
respondences. 

Lemma  6   For  any  function  v  G  H{div,K),  we  have: 

W4>  e  L\K),  I  ^V  -vdx^  I  <i)V  ■  vdx, 

Jk  Jk 

yd)  e  L^(dK),  I      4>n-vdsi=   [     4>n-vds:,. 

Jdk  JdK 
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Next,  we  state  another  Lemma  found  in  Raviart  and  Thomas  [30],  which  relates  the 
norms  and  semi-norms  of  scalar  and  vector  functions  and  their  correspondences. 


Lemma  7  For  any  integer  I  >  0,  we  have 

'{K} 


V</.  G  H\K),  \^\h'(K)  <  \\Bk\\^\Jk\  '^'|<^lH'(iC), 


Wie{H'{K))\  \v\^H>(K)y  <  \\BK\\'\\BK'\\\M^'\v\mK)y, 

where  \\Bk\\  and  \\BJ^^\\  denote  the  spectral  norm  of  the  matrices  Bk  and  BJ^  respec- 
tively. 

Construction  of  the  Raviart-Thomas  finite  element  spaces  using  the  spaces 
defined  on  the  reference  element  K. 

Let  r  be  a  non-negative  integer,  and  let  K  denote  the  reference  unit  triangle.  We 
outline  the  construction  of  the  Raviart-Thomas  spaces  Vh{^)  and  Q/i(f^),  of  order  r, 
on  f]. 

1.  Let  Vh,{K)  and  Qh{K)  be  finite  dimensional  linear  spaces  defined  on  the  unit 
reference  triangle  K . 

2.  Using  the  reference  spaces  Vh{K)  and  Qh{K).,  we  define  the  restriction  of  the 
Raviart-Thomas  spaces  Vh{^)  and  Qh{^)  to  any  non-degenerate  triangle  K  E 
r'',  which  we  denote  by  Vh{K)  and  Qh{K),  respectively.  We  define  Vh{K)  and 
Qh{K)  by  using  the  one-to-one  correspondences  obtained  by  using  the  affine 
linear  map  between  K  and  A',  as  discussed  in  the  previous  subsubsection.  i.e., 

Vh{K)  =  {v  e  H{div,K)  -.v^—^ie  Vh{K)}, 

Q,{K)  =  {4>e  L\K)  ■.4,^^4>^  Qk(K)}. 

3.  To  complete  our  construction  of  the  Raviart-Thomas  spaces  on  fi,  we  define  the 
spaces  by: 

F,(n)  =  {v^  H{div,  n)  :  ^Ik  €  V,{K)}, 

QHin)  =  {qeL\n):q\KeQH{K)}. 

In  the  next  subsubsection,  we  define  the  Raviart-Thomas  spaces  on  the  reference 
element  K. 
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Raviart-Thomas  spaces  Vh{K)  and  Qh{K)  defined  on  the  reference  element 
K. 

Here  we  define  the  Raviart-Thomas  spaces  Vh{K)  and  Qh{K).  We  remark  that  our 

notation  for  the  spaces  differs  from  that  used  in  Raviart  and  Thomas  [30],  with  the 

spaces  Vh{.)  and  Qh{-)  interchanged.  Let  r  be  a  non-negative  integer,  r  =  0, 1,2,3, 

The  space  Vh{K)  of  order  r,  is  a  finite  dimensional  linear  space  satisfying: 

(al)  {VriK)y  C  Vk{K); 

(a2)  dimension{Vh{K))  =  (r  +  l)(r  +  3); 

(a3)  ieVH{K)=^V-ieVr{K); 

(a4)  ieVh{K)=^K-J^Sr;    _ 

(a5)  Vf,{K)  n  H{div'>,  K)  C  {Vr{K)f; 

where 

•  Vt{K)  denotes  the  polynomials  of  degree  <  r,  on  K\ 

•  V  denotes  the  divergence  in  K\ 

•  h  denotes  the  unit  outer  normal  to  dK\ 

•  Sr  =  {s{x)  e  L\dK)  :  s{i)\i  e  Vr{e),  for  each  edge  e  C  dK}. 

The  following  Lemma  is  found  in  Raviart  and  Thomas  [30],  where  also  examples  of 
such  spaces  are  given. 

Lemma  8  For  a  space  Vh{K)  satisfying  conditions  (a2)  to  (o5),  a  function  v  G  Vh{K) 

is  uniquely  determined  by: 

(il)       the  values  of  v  •  n  at  {r  +  I)  distinct  points  on  each  edge  e,-  of  dK; 
(62)       the  moments  of  order  <  r  —  1  of  v,  i.e., 

We  define  the  Raviart-Thomas  pressure  space  Qh{K)  by: 

Qh{K)  =  Vr(K).  (1.17) 

Note  that,  using  the  definition  of  the  correspondence  for  scalar  functions,  we  obtain 
that  Qh{K)  =  Vr{K). 

REMARK.     Condition  (fcl)  can  equivalently  be  replaced  by  specifying  the  moments  of 
order  <  r  on  e,.  i.e., 

/   4){v  ■  n)dsi,    \I4>  e  Vr{ei), 

for  each  edge  e,  C  dK . 

19 


A  basis  for  Vh{K). 

Applying  Lemma  8,  we  construct  a  basis  for  V^(A')  as  follows.  First,  for  i  =  1,2,3 
and  j  =  0,...,r,  let  {</>,,;}  denote  an  orthonormal  piece- wise  polynomial  basis  for 
Sr{dK)  equipped  with  the  L^{dK)  inner  product.  On  each  edge  e^  C  dK,  let 
{^i.jYZo  ^  ^'■(<^i)  b^  t^^  Legendre  polynomials  of  degree  <  r  on  e,-  extended  by 
zero  on  the  other  edges,  e/,/  7^  i.  This  accounts  for  3(r  +  1)  basis  elements.  For  the 
remaining  r{r  +  1)  data,  we  use  the  nontrivial  distinct  moments 


/. 


k 
of  order  <  r  -  1  where  {pi]i^i       form  a  basis  for  {Vr-iiK))"^. 

To  construct  the  basis,  for  i  =  1,2,3  and  j  =  0,. . .  ,r,  we  let,  {uij}  denote  3(r  +  l) 
linearly  independent  functions  in  Vh{K)  satisfying: 

I  'yuVij  =  n  -Vij  =  4>x.j  on  dK, 

We  pick  the  remaining  r(r  +  l)  linearly  independent  basis  elements  {■u/}[i'^      in  Vh{K), 

as  follows: 

f  7„ti;  =  0,       V/, 

I  /a-  ^'  •  Pi^^  =  ^0'      "^^J  =  1,  ■  •  •  ,7-(r  +  1). 

Here,  S^^  denotes  the  Kronecker  delta  function, 

f     _  J    1      for  i  =  J 
^■^  -  \  0      for  2-  ^  j. 

Thus,  we  have  constructed  3(r  -f  1)  +  r{r  +  1)  linearly  independent  functions  which 
forms  a  basis  for  Vh{K).  We  will  denote  them  by  {i,,j}  U  {v-i},  for  i  =  1,2,3, 
j  =  l,...,(r  +  l),and/=  1, . . . ,  r(r  +  1). 

A  Dual  Basis  for  {Vh{K)y. 

Let  {Lij{.)}D{Li{.)}  denote  an  algebraic  dual  6asz5  with  respect  to  the  basis  {uij}  U  {iit} 
for  Vh{K),  i.e.,  the  {i,,j(.)}  U  {Li{.)}  denote  Unear  functionals  on  Vh{K)  C  C~(X), 
and  form  a  basis  for  the  dual  space  V/,(/f )'  satisfying: 

LiA^''.j')  =  ^iy^JJ'^  Vi,;,t',j', 

i^i,i(^/)  =  0,  Vt,j,Z, 

J^/(^.-,.)  =  0,  Vz,j,/ 

Li{ui-)  =  5i,/-,  V/,/'. 
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Thus  we  can  express  any  element  v  G  Vh{K)  by: 


».j 


REMARK .     We  can  express  the  dual  basis  explicitly  by, 

Je, 

and 


Li{v)  =  J^Pi-  ^<^^- 


An  example  of  the  lowest  order  Ravi  art -Thomas  spaces  defined  on  trian- 
gles. 

We  present  the  lowest  order  Raviart-Thomas  spaces,  i.e.,  for  r  =  0  on  triangles.  We 
have, 

Vh{K)  =  {(  ^;^^J^  ]  ■.ya,b,ceR;  wherex  =  (x:,i2)G^}, 

Qf,{K)  =  To{K). 
For  the  velocity  space  ^/.(Q),  we  have 

Vh{n)  =  {ve  H{div,n)  -.  v\k  =  {1/Jk)Bkv,  for  some  v  e  Vh{K)}. 

And  for  the  lowest  order  pressure  space  Qh{0.),  we  have 

QH{n)  =  {qeL'in):q\KeTo{K)}. 

Given  the  edges  of  the  triangles  forming  the  triangulation  t^  of  Q,  the  velocity  Uh  is 
uniquely  specified  by  its  normal  components  on  the  edges,  which  in  this  lowest  order 
case  is  constant  on  the  edges.  See  Figure  1.1.  A  basis  for  Qh{^)  is  easily  chosen;  the 
i"*  basis  function  is  one  on  the  i"'  triangle  and  zero  on  all  other  triangles. 

The  interpolation  map. 

We  end  this  subsection  by  introducing  an  interpolation  map  associated  with  the 
Raviart-Thomas  velocity  space  Vh{Q).  We  will  first  define  the  interpolation  map  fl 
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Figure  1.1:  Lowest  order  Raviart-Thomas  spaces  on  a  triangle. 


Uh  •  ^|ej  —  constant 


\/     Uh  •  n|e,  =  constcint 
e2\  '  ' 


Uh  ■  ^|ei  =  constant 


from  {H^{K))'^  onto  Vh{K).  We  discuss  briefly  why  this  is  a  continuous  map.  Then 
we  define  the  interpolation  map  U^  from  {H^{Q)y  onto  Vh{^),  by  letting 


h-r\ 


(nv 


J.'  < — » 


Uv,       \/K  G  r 


(1.18) 


We  need  to  check  that  this  is  a  valid  definition,  i.e.,  the  defined  function  is  in 
H{div,Q).  This  will  follow  by  construction,  from  the  compatible  normal  trace  cdong 
inter-element  boundaries,  see  Lemma  5.  The  interpolation  map  H^  will  also  be  shown 
to  be  continuous. 

DEFINITION  OF  fi.  By  Lemma  8,  given  a  vector  function  v  6  {E^{K)Y,  we  can 
determine  a  unique  interpolant  Uv  G  Vh{K)  satisfying: 

(cl)    V(^  e  Vr{e),      J^iUi  -  I)  •  ii4>dsi  =  0,     V  sides  e  C  dK, 
(c2)    V(p:,p2)€(7'._l(A'))^       4-(nt"-^).(pi,p2)i£  =  0. 

We  show  that  these  conditions  are  well  defined  for  {H^{K))^  functions  and  that 
the  interpolation  map  is  continuous.  We  use  the  property  that  for  functions  v  G 
{H\K)y,  its  normal  trace  -^J  G  H^'^-'{dK),  for  e  >  0.  This  follows  by: 

•  An  application  of  the  trace  Lemma  which  gives  77)  G  {H^^^{dK)Y. 

•  Using 

7„{;  =  n  •  7u. 

For  the  reference  triangle  K  the  unit  outer  normal  h  is  constant  on  each  edge 
e.  Thus  for  each  edge  e,  we  have 

7n^|e    G    H''\k). 


22 


•  Using  the  definition  of  H^^^  ^{^K)  on  the  one  dimensional  boundary  dK  to  get 
the  result. 

REMARK.  Note  that,  a  constant  vector,  which  is  a  function  in  {H^{K)y,  has  a  piece- 
wise  constant  normal  trace,  since  the  normal  is  discontinuous  across  the  vertices  of 
the  triangle  K.  Using  the  definition  of  the  E^'^{dK)  norm,  we  see  that  piecewise 
constant  functions  on  dK  are  not  in  H^'^[dK),  due  to  unbounded  norms. 
REMARK.  Note  that,  jn^^ulaK,  by  the  definition  of  El'',  is  the  affine  transformation 
of  a  polynomial  of  degree  less  than  or  equal  to  r,  on  each  edge  e  C  dK.  Also,  by  the 
definition  of  the  interpolation  map  the  inner  product  of  the  normal  trace  of  the  inter- 
polant  III"  with  functions  in  5,  equcJs  the  inner  product  of  the  same  functions  with 
the  normal  trace  of  u.  This,  together  with  Lemma  6  about  correspondences,  imply 
that  on  the  boundary  between  adjacent  elements  the  normal  trace  of  the  interpolant 
n''tr  is  compatible.  Thus,  the  interpolation  map  is  well  defined. 

Now  we  note  that  we  can  also  express  the  reference  interpolation  map  in  terms  of 
the  basis  and  dual  basis: 

Uv  =  Yi  Li_^{v)uij  +  YLi{v)ui. 

i.i  I 

By  using  the  fact  that  the  normal  trace  of  [H^(K.)Y  functions  are  in  H^^^~'{dK),  the 
linear  functionals  in  the  dual  basis  can  be  shown  to  be  continuous  on  {H^{K)y,  and 
thus,  n  is  continuous  on  {H"^ {K)Y .  This,  together  with  the  correspondences  between 
the  various  norms  involved  on  K  and  K  imply  the  continuity  of  II'',  cf.  Raviart  and 
Thomas  [30]. 

Lemma  9  Let  u  £  [H^{0.)Y ,  and  let  II''  denote  the  interpolation  operator  defined  by 
equation  1.18.   Then,  there  exists  a  positive  constant  c  such  that: 

We  present  details  of  the  proof  later  in  this  Chapter,  when  we  prove  the  continuity 
of  the  interpolation  map  on  H{div°,Q)  fl  {H^{Q)y. 

1.4.2      Quadrilateral  elements. 

In  this  section,  we  discuss  the  Raviart-Thomas  spaces  on  triangulations  consisting  of 
quadrilateral  elements.  See  Raviart  and  Thomas  [30]  and  Thomas  [32].  Let  K  denote 

23 


the  unit  reference  square  with  vertices, 

h,  =  (0,0);   <i2  =  (l,0);   a3  =  (l,l);   a4  =  (0,l). 

As  in  the  pre\'ious  section,  the  Ravi  art -Thorn  as  spaces  Vh{Cl)  and  Qhi^)  are  con- 
structed on  Q,  as  H{div,D,)  and  X^(Q)  subspaces,  respectively,  using  the  reference 
spaces  Vh(K)  and  Qh{K),  and  the  invertible  affine  linear  map  Fk  =  Ekx  +  Ik-  Since 
affine  linear  functions  map  squares  onto  parallelograms,  we  assume  that  the  elements 
K  in  the  triangulation  r^  of  Q,  are  parallelograms. 

Let  r  be  a  non-negative  integer.  We  present  the  definition  of  the  reference  spaces, 
Vh{K)  and  Qh{K)  of  order  r.  First,  for  non-negative  integers  7n,n,  let 

m       n 

VmJK)  =  {^  ^c„x;x^2  :  C'j  ^  ^;  where  i  =  (^1,12)  e  K}. 

We  let 

V^K)  =  {({>:, 1)2)  :  I'l  G  n  +  i,.(^),{;2  G  Vr,+i{K)}, 

and 
Note  that 

•  i.  e  v^iK)  =^  V  ■  i  G  g^(^)  =  n,.( A"). 

•  u  G  V'/,(A')  =>  7nUje,  G  ^r(e.),  for  each  edge  e,  C  OK. 
The  following  Lemma  is  stated  in  Raviart  and  Thomas  [30]. 

Lemma  10   A  function  u  G  Vh{K)  is  uniquely  determined  by: 

(1)  the  values  of  fnU  at  (r  +  1)  distinct  points  on  each  side  Ci  C  K , 

(2)  the  moments 

I  uix\x{dx,     0<r<r-l,     0<j<r, 

(  U2x\x\dx,     0  <r  <r,     0  <  j  <  r  -  I. 
J  k 

Condition  (1)  can  equivalently  be  replaced  by  the  moments  of  order  <  r  on  each  edge 

Ci  C  dK: 

f  hnudsi,    V<^G7'.(e,),    i=  1,2,3,4. 
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As  for  the  case  of  triangular  elements,  we  define 

Vh{K)  =  {v  :  V  •( »  V,  where  v  G  Vh{K)}, 

and 

Qf,{K)  =  {<l>:4>^^4>,  where  4>  £  Qh{K)}. 

Finally,  we  define 

V,{n)  =  {ve  Hidiv,n)  -.  €\k  e  Vh{K),VK  G  r''}, 
and 

Most  of  the  results  holding  for  the  triangular  case  holds  also  for  the  quadrilatercJ 
case. 

We  present  an  example  of  the  lowest  order  Raviart-Thomas  spaces  based  on  a 
rectangular  mesh. 

V^,{K)  =  {(  ^^^^'   J  ■.Va,b,c,deR;  where  x  =  (xj,!:)  e  ^^j, 

and 

Q^,{K)  =  Vo{K). 

As  stated  in  Lemma  10,  the  degrees  of  freedom  of  a  function  in  Vh{K)  can  be  specified 
by  the  values  of  {•  •  n  at  the  midpoint  on  each  of  the  four  sides  of  K.  Of  course,  v  •  h 
is  constant  on  each  edge,  for  the  lowest  order  case.  See  Figure  1.2. 

1.4.3      Convergence  results  and  other  properties. 

We  list  some  of  the  properties  of  the  Raviart-Thomas  spaces  Vh{Q)  and  Qh{^),  which 
we  use  later  on. 

The  functions  in  V/i(Q),  which  are  discrete  divergence  free,  are  divergence  free  in  the 
L^{Cl)  sense,  i.e.,  if  Vh  G  V/,(fi)  and 

/  <^(V  •  Vh)dx  =  0,      V<^  €  Qh{Cl),  then  V  ■Vh  =  0. 
Jn 

PROOF.     From  the  definition  of  the  reference  spaces,  we  see  that, 

Vf,eVh{K)^V-ieQH{K). 
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Uh  -Tile, 


64 

constant 


Uh  ■  ■^lej  =  constant 


ea 


PkIk  —  constant 


ei 


62 


Uh  •  n L    =  constant 


Uh  ■  7i|ej  =  constant 


Figure  1.2:  Lowest  order  Raviart-Thomas  spaces  on  rectangles. 


Thus, 

/  <^(v  •  4)^2;  =  0,    y(f>eQh{n) 

Ju 
=^   f  ct>(V-  VH)dx  =  0,       V<^  e  Qh{K) 

JK 

=^.  [  4>V.  indi  =  0,       V</.  G  ^/.(r). 
By  substituting  ^  =  V  ■  Vh,  we  obtain  V  ■  i'h  =  0. 

i.e.,  (V  •  Vh)  —  0.  Hence,  for  the  Raviart-Thomas  spaces,  a  discrete  divergence  free 
function  is  divergence  free  in  L'^{Q,).  □ 
An  additional  property  is: 

If  n'^i^  is  defined  for  a  H{div°,n)  function  v,  then  U^v  6  H{div°,n).  i.e.,  the  in- 
terpolation map,  when  defined,  takes  divergence  free  functions  into  divergence  free 
functions. 

PROOF.  See  Raviart  and  Thomas  [30].  We  show  that  when  the  interpolation  map 
is  well  defined  for  a  divergence  free  function  u  its  interpolant  U^u  is  also  divergence 
free.  Consider 

/  qV  ■  Uvdx  =   Iv-  nqdsi  -  f  Uv  ■  Vqdx.  (1.19) 

JJK  JdK  JK 

For  q  G  Vt{K),  for  triangles,  (or  q  £  Vr,r{K)  for  the  quadrilatercJ  case,)  we  obtain 

that  Vq  e  {Tr-\(K)f,  (and  Vg  G  Vr-i,r(K)  x  Vr,r--L(K),  for  the  quadrilateral  case.) 

For  the  same  q,  we  obtain  that  its  restriction  to  any  edge  e  in  the  boundary  is  a 

polynomizJ  in  'Pr(^).  Both  these  imply  that  we  can  replace  fiu  by  u  in  equation  (1.19), 
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by  using  the  definition  of  the  reference  interpolation  map  fl.  By  an  application  of 
the  divergence  theorem,  we  obtain  that  this  expression  equals  zero,  since  V  •  u  =  0. 
Thus,  the  interpolant  is  discrete  divergence  free  which  implies  divergence  free  in  the 
L^  sense. 

We  now  check  whether  the  appropriate  discrete  subspaces  for  the  Neumann  bound- 
ary value  problem  (1.7)  satisfies  the  ellipticity  and  uniform  inf-sup  conditions.  Of 
course,  this  would  be  a  sufficient  condition  for  the  convergence  of  the  discrete  prob- 
lem. REMARK .  For  Neumann  boundary  conditions,  the  appropriate  discrete  subspaces 
of  Ho,dn{div,Cl)  and  L^{Q)/R  used  in  discretising  problem  (1.7)  cire 

x,(n)  =  vj.(n)  n  HoMd^^'^^)^  and  n(n)  =  Q^,{n)  n  (i'(fi)/iZ), 

respectively,  where  V/,(f2)  and  Qhi^)  denote  the  Raviart-Thomas  finite  element  spaces 

of  order  r  based  on  the  triangulation  r  . 

REMARK.     Note  that  since  a(.,.)  is  Ho,dn{(iiiJ°,n)- elliptic  with  a  constant  a,  so  is  its 

subspace  Xh{^),  with  the  same  constant  a.  See  equation  (1.13). 

REMARK.     A  lemma  in  Raviart  and  Thomas  [30]  states  that  given  q^  G  Qh,  there 

exists  Vh  E  Vh{Cl)  such  that 

^■"^^     "    ^'"  (120) 

where  c  >  0  is  a  constant  independent  of  h.  This  is  equivaJent  to  a  uniform  inf-sup 
condition  for  Dirichlet  boundary  value  problems.  However,  for  Neumann  boundary 
value  problems,  given  qn  G  Qhi^)  H  L^{Q.)/R,  we  need  to  find  Vh  satisfying  the  above 
bounds,  but  with  4  G  Xhi^)  =  Vh{i})  H  Ho_aQ{div,n).  See  Thomas  [32].  We  give  a 
proof  using  an  extension  theorem  for  the  Raviart-Thomas  velocity  space  V/i(fi).  See 
Theorem  2,  discussed  in  the  next  section. 

By  equation  (1.20),  given  qn  E  Qh{^)f^L^{^)/R,  there  exists  Vh  G  F/,(fi)  satisfying 
the  above  bounds.  But  jn't^h  ^  0,  in  general.  However,  since  V  -Vh  =  9^,  we  have  by 
an  application  of  the  divergence  theorem  that  Jg^  -ynV^dsx  =  0.  Now,  by  an  extension 
theorem  discussed  in  the  next  section,  there  exists  E^jn^h  G  Vh{^)r\H{div°yQ)  such 
that  V  •  (E^fr^Vh)  =  0  and  such  that 


\E^'rnVh\\n(d^v.n)  <  c{^)\\'rr.Vh\\H-u^any 
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Here  c{Q.)  is  a  positive  constant  independent  of  h.  By  an  application  of  the  normal 
trace  lemma,  we  obtain 

Thus  letting  {1^  =  v^  -  E''-fnVh,  we  obtain 

V-TXh     =     qh, 

where  c  >  0  is  another  constant  independent  of  h.  Note  that,  fnUh  =  0.  This  shows 
that  the  divergence  map 

is  an  isomorphism  from  Vh{n)  n  HoMdiv,n)  onto  Qh{n)  n  i2(Q)/-po(n).  This  is 
equivalent  to  the  uniform  inj-sup  condition^  by  one  of  the  remarks  following  Theo- 
rem 1. 

Approximation  error  for  the  Raviart-Thomas  finite  element  spaces. 

The  following  result  is  found  in  Raviart  and  Thomas  [30]. 

Lemma  11  Suppose  that  p  e  H'''^^{Q.)  and  Ap  G  i/'"^^(Q),  for  some  integer  r  >  0. 
Then  the  discretisation  error  using  the  Raviart-Thomas  finite  element  spaces  of  order 
r  IS 

<  ch""^  (|p|H'+'(n)  +  |pl//'+2(n)  +  |Ap|H'+"(n)), 
for  a  positive  constant  c  independent  of  h. 

This  result  implies  convergence  of  the  mixed  finite  element  method  using  the  Raviart- 
Thomas  finite  element  spaces  when  the  continuous  problem  has  sufficiently  smooth 
solutions. 
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1.5      An  extension  theorem  for  Raviart- Thomas  fi- 
nite element  spaces. 

In  this  section,  we  introduce  an  extension  theorem  for  the  Raviart-Thomas  velocity 
spaces,  which  will  be  used  later  to  estimate  the  rate  of  convergence  of  various  domain 
decomposition  algorithms.  By  an  extension  map,  we  mean  a  continuous  map  from 
a  space  of  boundary  values  on  dD.  to  a  related  space  of  functions  defined  on  Q. 
As  before,  we  let  Vh{Cl)  denote  the  Raviart-Thomas  velocity  space  of  order  r  on  a 
uniformly  shape  regular  triangulation  r^  of  Q.  Note  that  the  normal  trace  of  elements 
in  Vh{Cl)  need  not  be  zero.  Let  us  denote  by  Vh{dQ,),  the  space  of  normal  traces  of 
elements  in  Vh{^)  : 

Vh{^n)  =  {fr.v■.VHev^,{n)}. 

Thus 

and  by  applying  the  normal  trace  lemma  to  v^  £  Vh{^)  we  obtain: 

||7„4||i/->/^(an)  <  Ci^)\Wh\\H(d,v,Q)- 

We  first  note  that  since  7^  is  continuous  and  onto  Vh{dn),  by  an  application  of  the 
open  mapping  theorem,  7^  has  a  continuous  right  inverse.  This  right  inverse  provides 
an  extension  of  the  normal  trace  boundary  values  in  Vh{dCl)  into  the  Raviart-Thomas 
velocity  space  Vh,{Q).  However  the  norm  of  this  right  inverse  may  depend  on  h. 

When  we  apply  extension  theorems  to  domain  decomposition  methods  to  estimate 
the  rate  of  convergence  of  these  iterative  methods,  we  obtain  theoreticcJ  bounds  which 
are  proportional  to  the  norm  of  the  extension  map.  Thus,  if  for  each  h,  we  could  find 
an  extension  E^,  whose  norm  is  bounded  independent  of  h,  then  our  bound  for  the 
rate  of  convergence  of  the  iterative  method  would  be  independent  of  h.  In  another 
application,  a  proof  of  the  uniform  inf-sup  condition  for  Neumann  boundary  value 
problems  for  mixed  formulations  of  elliptic  problems  using  Raviart-Thomas  elements, 
also  requires  an  extension  whose  norm  is  bounded  independent  of  h. 

For  ^^-conforming  finite  element  functions,  such  extensions  or  prolongations  from 
the  finite  element  trace  spaces  to  the  finite  element  spaces,  uniformly  bounded  in  /i, 
have  been  constructed  in  several  ways.  See  Astrakhantsev  [1],  Widlund  [34],  Bj0rstad 
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and  Widlund  [4],  Bramble,  Pasciak  and  Schatz  [5],  and  Matsokin  [25].  Using  ideas 
from  Bj0rstad  and  Widlund  [4],  and  Bramble,  Pasciak  and  Schatz  [5],  we  construct 
a  uniformly  bounded  extension  for  H{div°,Q,)  conforming  Raviart-Thomas  velocity 
spaces. 

We  introduce  some  notations.  Let 

Sh{n)     =  VH{n)nH{div°,ri), 
Sh{dn)  =  {jnVH-.VHeSHin)}. 

Note,  that  equivaJently,  we  have 

ShiOn)  =  {g,  G  VnidQ)  :  /    g^ds,  =  0}. 

Jan 

Before  we  describe  the  extension  theorem,  we  state  a  few  Lemmas  that  will  be 

used  in  its  proof.  Let  K  denote  the  reference  element  and  e,-  C  OK  one  of  its  edges. 

Lemma  12  LdO  <  e  <  1,  and {1  ^  {H'{K))^r]H{div°,K).  Then,    7„u|e-.  G  H'-^^\ei). 
Furthermore 

\\lnU\\H'-^'^e.)   <   C(e,iir)||ti||(^.(^.)),, 

for  some  positive  constant  c{e,K). 

PROOF  OF  LEMMA.  By  an  application  of  the  trace  lemma,  since  the  normal  n  on  the 
edge  e,-  is  constant,  we  obtain  that 

is  a  continuous  map.  Similarly,  by  an  application  of  the  normal  trace  lemma,  and  the 
fact  that  restrictions  of  H~^^'^{dK)  functions  to  edges  e,  C  OK  are  in  (.ffoo  {^i))'i  '^^ 
obtain  that 

^^■.H{dzv\K)—.{Hll\iO)\ 

is  a  continuous  map.  See  the  remark  in  Chapter  4,  section  2,  for  more  details  about 
this.  Thus,  we  have  that 

fr.:{L'{K)ynH{div°,K)      —.     {Hll\e,)y, 
fr.:{H\K)ynH{div\K)    -^    H'/'{e,), 

are  both  continuous  maps.  Using  the  theory  of  interpolation  spaces,  we  obtain  that 
the  following  map  is  continuous: 

'^r.:{H'(K)fnH{div\K)—.H'-'l\e,)       for  0  <  £  <  L      •'' 
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i.e.,  there  exists  a  positive  constant  c{e,K)  such  that 

ll7ni||H-'/2(e.)  <  c(e,iir)||'u||(^.(^))j.D 

Lemma  13  Let  e  >  0  and  let  Cl  be  a  2  dimensional  polygonal  region  triangulated  by  a 
family  o/uniformly  shape  regular  triangulations  {r^}.  Let  u  6  {H'{Cl)yr\H{div°,Q,), 
and  let  U.^  denote  the  interpolation  operator  defined  onto  the  Raviart- Thomas  finite 
element  velocity  space  Vh{Q.).  Since  u  is  divergence  free,  n'*Tr  G  5^(0)  =  V/j(f2)  fl 
H{div°,Q,),  and  there  exists  a  positive  constant  c(e,Q)  such  that: 

PROOF  OF  LEMMA.  See  the  section  on  Raviart-Thomas  elements  for  the  definition 
of  the  interpolation  map.  We  first  consider  the  interpolation  map  on  the  reference 
element  K.  By  using  the  representation  of  the  interpolation  map  as  a  sum  involving 
linear  functionals  and  its  basis,  we  first  check  that  the  linear  functionals  forming  the 
dual  basis  are  well  defined  and  continuous  for  {H^{K)Y  Pi  H{div°,K)  functions.  By 
Lemma  12,  7n^'|e.  S  H'~'^^'^{e,).  Thus,  since 

Liji'^')  =  /    4>i.j'^  -vdsi, 

Je, 

and  since  4>i.j  ^  ^""(e,),  we  obtain  the  bound 

\Li^j{v)\  <  c(e,A')||7„i)||H.-i/2(.-.)||<^i,j|lH>/=-'(e.)   <  c(e,-ft')||i'i|(H«{K))sll<^.jllH'/3-.(e.)- 

Similarly,  since 

Li{v)  =     .  Pi  •  vdx, 

J  K 

and  since  pj  G  (T'r-iC-ft'))^,  ^^  ^^"^^ 

\Li{v)\  <  \m{LHK))^\\P'\\iL^K)y  ^  \M{H'{K)r\\pt\\{L^f:)y 
By  combining  these,  we  obtain, 


»,j 
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Here  c{e,K)  depends  on  the  sum  of  the  bounds  for  each  linear  functional.  Thus,  fl 
is  a  bounded  map. 

Next,  we  discuss  why  II''  is  a  well  defined  and  bounded  map.  Recall  that,  on  each 
element  K, 

where  v  * >  v\k-  We  recall  that,  by  Lemma  6 


Thus,  since  v  e  H{div,Q),  and  since  the  interpolant  7„(n''i7)|e  6  P,_i(e),  we  deduce 
that  U'^v  G  H{div,n).  To  show  that  the  interpolation  map  II''  is  bounded,  we  observe 
that: 

1.  Using  the  transformations  of  semi-norms  and  norms  between  K  and  K,  as 
given  in  Lemma  8,  and  the  boundedness  of  the  reference  interpolation  map  H, 
we  obtain: 

since,  \\Bk\\  =  0(/ia-);  |Ja'|  =  0{h^j^),  etc. 

2.  Summing  up  the  above  result  over  VA'  €  t'',  we  obtain  that  ||n'*||  <  oo.  □ 

Lemma  14  Let  e  >  0  and  let  ^  be  a  2  dimensional  polygonal  region  triangulated  by  a 
family  o/uniformly  shape  regular  triangulations  {t^}.  Let  u  G  {H'{Q)Yr\H{div°,Q,), 
and  let  II''  denote  the  interpolation  operator  defined  onto  the  Raviart- Thomas  finite 
element  velocity  space  Vh{Cl).    Then,  there  exists  a  positive  constant  c(€,  Q)  such  that: 

W  -  n')^b(,.„o,n)  <  c{e,n)h^\\u\\^jj.m^. 

Here  I  denotes  the  identity. 

PROOF  OF  LEMMA .  For  e  >  0,  by  Lemma  13,  we  have  that  II''  is  a  bounded  linear  map 

n' :  {H'{n)f  n  H{div'',n)  -^  s^in). 
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Using  this  and  a  quotient  space  scaling  argument,  we  will  show  that 

\\{I-U')uy^,,^o,^)<c{e)h'\\u\\^H.^a)y. 

We  now  present  the  details  of  this  quotient  space  argument,  which  is  also  used  in  the 
study  of  the  convergence  of  the  finite  element  method.  See  Ciarlet  [12],  and  Girault 
and  Ravi  art  [15]. 

Let  K  denote  the  reference  unit  element,  a  triangle  or  a  square  as  the  case  may 
be,  and  let  K  be  an  element  in  the  triangulation  t^  of  Q.  Let  Fjc  denote  the  adfine 
linear  map  of  K  onto  K. 

FK-.ieK  —^  Fk{x)  =  Bkx  +  fcA', 

where  Bk  is  a  2  x  2  invertible  matrix  and  h^  is  a  2- vector.  Let  Jk  =  det[BK)-  We 
use 

u  < >  u 

to  denote  the  correspondences  between  vector  functions  on  K  and  on  K,  as  described 
in  an  eariier  section. 

By  the  relation  between  corresponding  vector  functions  on  K  and  K ,  we  obtain: 

||(7-n'')tiii(i,(K)p  <  |Jicr^/'||B;,||  ||(/-n)5||(^,(^)),. 

For  e  >  0,  by  the  boundedness  of  I  and  E,  we  obtain: 

\JKr"\\BK\\    \\{I  -"^M^mkW  <\Jk\-'^'\\Bk\\    11/  -fill    H(H.(X))- 

Since  both  the  identity  /  and  the  interpolation  map  fl  preserve  constants,  i.e.  elements 
in  (Po(^))^  we  may  replace  \\u\\(^ij.^k)Y  ^y  the  quotient  space  {H'(K)f  l{Vo(K)f 
norm.  By  Poincare's  inequality  this  is  bounded  by  the  {H^{K)y  semi-norm.  Thus, 
we  have 

||(7-n'')tr||(,.(K))^  <c(^)1Ja|-^/^||B;,||  ||7-n||M(H.(^„,. 

Now,  mapping  u  back  to  u  on  K  and  using  bounds  on  the  Jacobian,  we  get. 
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Since  Cl  is  assumed  to  be  a  polygonal  domain,  we  may  apply  the  elliptic  regularity 
theory  for  this  Neumann  problem  only  for  0  <  e  <  1/2,  see  Necas  [27],  Grisvard  [19], 
and  Babuska  and  Aziz  [3].  If  c//>  G  H"^l'^{dQ.)  then  ip  G  E^+'{^).  Furthermore  we 
have 

llvllH'+Mn)  ^  c(n,£)||g/,||H-'/2(an)5 

for  some  positive  constant  c{Q.,(.).  Note,  gh  G  L^{dD,),  therefore  gh  G  H*~^^^{d^) 
since  the  normal  trace  space  Sh{dQ)  consists  of  piece- wise  polynomial  functions  with 
discontinuities  at  the  vertices.  If  Q,  were  a  smooth  domain,  the  elliptic  regularity 
theory  would  hold  for  all  e  >  0.  For  the  special  case  of  e  =  0  these  results  can  also  be 
obtained  by  using  the  standard  H^{Cl)  variational  formulation.  For  details  see  Necas 
[27],  Grisvard  [19],  and  Babuska  and  Aziz  [3]. 

Because  Eg^  =  Vy?,  it  satisfies  V  •  Eg^  =  0.  Therefore  Eg^  G  H{div°,Q,).  We  also 
have 

\\Egh\\{n'(Q]y  <  c(fi,e)||5h||H«-'/2(sn)5 

for  some  positive  constant  c{Q,e). 

However,  Eg^  need  not  be  in  5^(0).  In  order  to  obtain  an  extension  E'^gn  in  the 
finite  element  space,  we  use  the  interpolation  map  onto  the  Raviart-Thomas  velocity 
space,  i.e.,  we  let 

E^'gH^U^Eg,), 

where  II''  denotes  the  Raviart-Thomas  velocity  interpolation  map.  See  section  1.4, 
on  the  Raviart-Thomas  finite  element  spaces  for  the  definition  and  properties  of  this 
interpolation  map.  Now,  we  use  Lemma  14,  which  states  that: 

Since 

U'Eg,  =  Eg,  +  {U''Eg,-Egh), 
by  an  application  of  the  triangle  inequality,  we  have 

P''^5/.|l5(di„o,n)  <  \\Egh\\i}^d,^o^Q)  +  c{€,n))h'\\Egh\\(n'(n)y-  ■ 
Using  the  elliptic  regularity  results  we  have 
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By  using  the  equivalence  of  all  norms  on  a  finite  dimensional  space,  we  can  bound 
the  H'~'^^^{dQ,)  norm  of  gh  in  terms  of  its  H~^^^{dO,)  norm,  with  a  constant  possibly 
depending  on  h.  We  estimate  this  dependence  on  h,  by  using  an  inverse  inequality 
for  the  finite  element  space  Vh{dQ,).  We  obtain  that 

\\9h\\B'-'/\da)  <  c(£,f2)/i"'||5/,||5-i/2(an).  (1.21) 

Though  Vh{dQ)  consists  of  piece-wise  polynomial  functions  which  are  discontinuous 
across  vertex  nodes,  we  can  still  prove  such  an  inverse  inequality.  We  outline  how 
this  can  be  shown.  First,  we  prove  inequality  (1.21)  for  small  e  >  0  with  —1/2  re- 
placed by  -1-1/2.  Following  that,  we  use  duality  and  interpolation  theory  to  obtain 
intquality  (1.21).  For  details  about  this,  see  Babuska-and  Aziz  [3].  To  prove  inequal- 
ity (1.21)  with  —1/2  replaced  by  1/2,  we  use  formula  (1.2)  to  express  the  norms. 
Then,  we  partition  dQ  into  a  disjoint  union  of  edges,  each  of  which  is  the  the  edge 
of  some  element  in  the  triangulation  r'*.  Thus,  the  norms  can  be  decomposed  as  as  a 
sum  of  terms  involving  the  edges  that  constitute  dQ.  We  then  consider  the  quotient 
of  the  semi-norms;  this  can  be  estimated  by  mapping  the  numerator  and  denominator 
onto  a  reference  unit  element  by  using  an  affine  linear  map.  Since  gh  restricted  to 
any  edge  is  a  polynomial,  we  can  bound  the  quotient  on  the  reference  element  by  a 
constant  depending  only  on  the  dimension  of  the  polynomial  space  of  the  reference 
edge,  and  the  quotient  of  semi-norms  on  this  polynomial  space.  Using  the  bounds  for 
the  Jacobian  of  the  map  between  the  edge  and  a  reference  element,  in  terms  of  h,  we 
can  establish  inequality  (1.21)  with  —1/2  replaced  by  +1/2.  We  do  not  present  the 
details.  See  Babuika  and  Aziz  [3],  Ciarlet  [12]  and  Girault  and  Raviart  [15].  Putting 
these  resiilts  together,  we  obtain 

P''-E£?/»||^(<j,„o,n)  <  c\\gh\\H--^/^aQ)  +  c(e,f2)/i'/i~'||5,,||^-i/j(an). 
Thus,  £"'*  =  U^E  is  bounded  independently  of  h 

Note  that  by  construction,  II''£'^;;  G  5/,(f2),  and  -f^Egh  =  gh-  Thus,  E^  =  U'^E  is 
a  uniformly  bounded  extension.  □ 
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1.6      Iterative  methods  for  solving  symmetric,  pos- 
itive definite  linear  systems. 

Let  A  be  an  real,  nxn  symmetric,  positive  definite  matrix,  and  b  6  i2".  We  consider 
solving  linear  systems  of  the  following  form 

Find  X  e  R",  such  that       Ax  =  b, 

using  iterative  methods. 

Let  M  be  a  symmetric,  positive  definite  matrix,  which  is  easier  to  invert  than 
A,  in  the  sense  of  solving  Mx  =  b.  For  i  =  0, 1,2, . . .,  given  the  iterate  i',  the  next 
iterate,  x'"*"^  is  defined  by 

Mx'+^  ={M  -  A)x'  +  b.  (1.22) 

The  initial  iterate,  i°  is  arbitrary  and  can  be  chosen  to  be  zero.  We  note  that 
equation  (1.22),  when  considered  as  a  mapping  from  x'  to  x'"'"\  has  a  fixed  point  at 
I,  i.e.,  X  is  the  unique  solution  of 

Mx  =  {M  -  A)x  +  b.  (1.23) 

Subtracting  equation  (1.22)  from  equation  (1.23),  we  obtain, 

M{x  -  x'+')  =  {M  -  A){x  -  x^). 

Letting  e'  =  x  -  x\  denote  the  error  in  the  iterate  x\  we  obtain  Af  e'"^^  =  {M  —  A)€\ 
i.e., 

e'+^  ={I  -  M-'A)e\ 

Let  p{I  —  M~M)  denote  the  spectral  radius  of  7  —  M~M.  Using  the  eigenfunction 
expansion  of  e°,  we  see  that  if  p{I  —  M~M)  <  1,  the  iterative  method  converges.  We 
note  that  I  —  M~^A  is  diagonalisable  since  it  is  similar  to  M^^^{I  —  M~^A)M~^^^  = 
{I-M~'^^^AM~'^^^),  which  is  a  symmetric  matrix.  In  the  rest  of  the  thesis,  we  consider 
iterative  methods  to  solve  the  saddle  point  problems  obtained  from  the  mixed  finite 
element  discretisation  of  elliptic  problems  using  the  Raviart-Thomas  spaces. 

For  various  choices  of  M,  the  spectral  radius  p{I  -  M~^A)  may  not  be  less  than 
one.  In  such  cases,  we  can  use  an  "acclerated"  iterative  method,  such  as  the  conjugate 
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gradient  method,  to  obtain  convergence.  We  may  also  use  the  acclerated  methods  even 
when  the  basic  iterative  method  is  convergent.  We  breifly  describe  the  preconditioned 
conjugate  gradient  algorithm,  and  then  present  the  error  propagation  results.  For 
details,  see  Axelsson  and  Barker  [2],  Golub  and  Van  Loan  [17],  and  Hageman  and 
Young  [20].  The  algorithm  is  described  for  the  symmetric,  positive  definite  linear 
system  Ax  =  b,  with  a  symmetric,  positive  definite  preconditioner  AI .  Let  e  >  0 
denote  a  tolerance  used  as  a  stopping  criterion,  and  x  the  iterate  obtained  when  the 
stopping  criterion  is  satisfied. 


x°  =  0 

r°  =  6 

Forz'=  1,2,... 

if  ||r'-i  1  <  |jr°||e  then 

Set  X  =  x'"''  and  quit. 

else 

Solve  Mr'-^  =  r'-^  for  z^-\ 

^'  =  [2-^r-l]/[z'-^^'-2]     (/?' 

^0) 

p,  =  ,.-i^^y-:     (pi=,o^ 

a'  =  [z^'\r^-']/{p\Ap^] 

x'  =  x'-^  +  qV 

r'  z^r'-^  -  q'-V 

en 


dfor. 


Here  {x,y]  denotes  the  euclidean  inner  product  x'^y,  when  A  and  M  are  symmetric 
in  the  usual  sense.  If  M  =  /,  and  A  is  symmetric  in  an  inner  product  which  is  not 
the  euclidean  inner  product,  then  [.,.]  will  denote  that  inner  product. 

To  obtain  a  computationally  more  efficient  algorithm,  we  assume  that  it  is  easier  to 
invert  the  preconditioner  M  than  the  matrix  A,  in  the  sense  of  solving  linear  systems; 
i.e.,  it  is  computationally  easier  to  solve  Mz  =  d,  than  to  solve  Az  =  d.  Furthermore, 
we  assume  that  the  spectral  condition  number  of  M~^^^AM~^^^  is  smaller  than  the 
spectral  condition  number  of  A.  The  spectral  condition  number  of  a  nonsingular  matrix 
B  is  defined  to  be  ||B|]||5"^||,  where  ||B||  and  |1B~^||  denote  the  euclidean  norm  of  B 
and  B~^  respectively.  Note  that  the  algorithm  only  involves  use  of  the  product  Ax, 
and  the  solution  of  My  =  z,  and  other  vector  and  scalar  operations.  In  particular,  we 
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do  not  need  to  assemble  the  matrix  A  or  M,  we  only  need  to  know  how  to  compute 
Ax  and  M~'^z. 

The  error  for  the  conjugate  gradient  method  with  M  =  /  is  bounded  by 


fel)'li^-^°ll- 


|i  —  a^'IU  £  2  I  — ;= — -  I    III  —  X 


where  the  >l-norm  \\w\\a  =  (w^AwY^^,  and  k  —  k,{A)  is  the  spectral  condition  number 
of  A.  A  similar  result  holds  for  the  preconditioned  conjugate  gradient  method  with 
k{A)  replaced  by  k{M~'^^^AM~'^^'^).  Note  that  the  spectral  condition  number  of  j4  is 
equivalently  defined  by  A2/A1,  where 

x'^Ax 
Ai  <  -^^  <  A2, 

are  the  maximum  and  minimum  eigenvalues  of  A  given  by  the  extremum  of  the 
Rayleigh  quotient  of  A.  Similarly,  the  spectral  condition  number  of  M~^^^AM~^^^  is 
equivalently  given  by  Ao/Ai,  the  quotient  of  the  extremum  of  the  generalised  Rayleigh 

que  .lent 

x'^Ax 

x^Mx 
Finally,  we  remark  that  if  M  =  /,  and  if  A  is  symmetric,  positive  definite  in  an 
inner  product  [., .]  which  may  be  difi"erent  from  the  Euclidean  inner  product,  then  the 
same  conjugate  gradient  algorithm  holds,  with  use  of  the  [.,.]  inner  product. 
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Chapter  2 

The  Schwarz  methods  for 
Raviart-Thomas  elements. 


In  this  chapter,  we  introduce  the  multiplicative  and  additive  Schwarz  methods  which 
are  iterative  methods  to  solve  symmetric  positive  definite  problems.  First,  we  describe 
these  methods  in  a  Hilbert  space.  Following  that,  we  outline  their  application  to 
saddle  point  problems.  We  then  study  in  detail  the  application  of  these  methods  to 
the  particular  case  of  mixed  finite  element  discretisations  of  elliptic  problems  with 
Neumann  boundary  conditions  using  the  Raviart-Thomas  spaces.  Finally,  we  present 
the  results  of  various  numerical  experiments  conducted  using  these  methods. 

2.1      The  Schwarz  alternating  methods  on  a  Hilbert 
space. 

Let  T~i  he  &  finite  dimensional  Hilbert  space  with  inner  product  a(., .),  and  let  /(.)  be 
a  continuous  linear  functional  on  Ti.  Consider  the  following  problem: 

{Find  u  ^  7i  such  that  ,.-  ,. 

a{u,v)  =  f{v),  yv^n.  ^^^ 

In  matrix  language  this  is  just  a  symmetric,  positive  definite  linear  system,  and 
it  could  be  solved  by  a  standard  iterative  method  such  as  the  conjugate  gradient 
method.  We  first  describe  the  Schwarz  alternating  method,  or  as  we  wiU  refer  to  it, 
the  multiplicative  Schwarz  method. 
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2.1.1      The  multiplicative  Schwarz  method. 

The  multiplicative  Schwarz  method  is  an  iterative  method  to  solve  problem  (2.1),  cf. 
Lioii^  [22].  We  assume  that  the  space  7i  is  decomposed  as  a  sum  of  smaller  subspaces. 
Then,  during  each  iteration,  which  usually  consists  of  a  fevp  fractional  steps,  each  of 
which  involves  a  defect  correction  on  one  of  the  subspaces,  we  update  the  old  iterate 
so  that  the  iterate  in  the  current  fractional  step  has  zero  residual  with  respect  to 
the  appropriate  subspace.  Thus,  each  iteration  involves  the  solution  of  problems  of  a 
smaller  size  and  often  some  subproblems  can  be  solved  concurrently,  i.e.,  in  parallel. 
We  now  describe  this  in  more  detail.  For  i  =  0, . . .  ,n,  let  Tii  be  subspaces  of  7i 
that  sum  to  H: 

Let  P,  denote  the  orthogonal  projection  onto  7ii,  in  the  a(.,.)  inner  product.   Thus 
PiU  G  Tii  satisfies 

a{P,u,v)  =  a{u,v)     Vv  G  7i,.  (2.2) 

Finding  P,k  is  equivalent  to  solving  a  "subproblem"  on  7i,. 

Let  u°  G  "W  be  an  arbitrarily  chosen  initial  iterate.  Given  the  i"*  iterate,  ti',  for 
i  =  0,1,2,...,  we  compute  the  next  iterate,  u*"'"^,  in  n  +  1  fractional  steps.  Each 
fractional  step  consists  of  a  defect  correction  on  one  of  the  subspaces.    i.e.,  for  j  = 

-1.2+1 

0. . . .  ,n,  let  u'    "+'   G  "H  be  sequentially  defined  by: 


where  6u     "+'  G  Tlj  is  chosen  so  that 


This  step  is  known  as  the  defect  correction  on  Tij.  Thus,  ^u'"*""*'  =  Pj{u  —  ■u'''""+'), 
is  the  solution  of 

a{Su'-^^,v)  =  a{u-u'-^^,v)  =  f{v)-a{u'-^^,v)    VveTij.  (2.3) 

Since  the  subproblems  are  of  a  smaller  size,  they  may  be  solved  with  more  ease 
than  the  full  problem.  In  case  some  subspaces  {7ij}  are  mutually  a(.,  .)-orthogonal, 
then  we  can  solve  some  of  these  subproblems  concurrently,  i.e.,  in  parallel.  We  will 
discuss  this  in  the  applications.  Note  that  the  right  hand  side  in  equation  (2.3)  can 
be  computed  without  knowing  u,  since  a[u, .)  =  /(.)  is  given. 
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The  error,  e'  =  u  —  n\  is  found  to  satisfy 

p'"'"^  —  P-"-  P-"-     . . .  P-^  P-p' 
^        —  -r i  -r^fc-i  •  •  •  M  -To  e  , 

where  Py^  =  I  —  P^  denotes  the  a(.,  .)-orthogonal  projection  onto  Tif .  Thus,  we  have 

lle'li  <  p'\\e\       where  p  =  \\P,^P,t,  •  •  •  P,^P,^ ]|, 

and  the  norm  is  defined  by 

||iy||  =  s/a{w,w). 

Since  orthogonal  projections  have  norms  bounded  by  one,  p  <  1.  However,  for 
convergence  of  this  iterative  method,  we  need  p  <  1. 

The  following  Lemmas,  discuss  sufficient  conditions  for  the  convergence  of  the 
multiplicative  Schwarz  method.  The  main  condition  is  that  the  subspaces  Tij  sum  to 
Ti.  However,  explicit  bounds  for  the  convergence  factor,  p,  can  be  obtained  in  terms 
of  other  parameters. 

The  next  Lemma  establishes  certain  inequalities  satisfied  by  the  projections.  It  is 
also  used  in  the  proof  of  convergence  of  the  additive  Schwarz  method.  It  is  given  in 
Lions  [22]  for  the  case  k  =  2.  Also  see  Dryja  and  Widlund  [131. 

Lemma  15  (Lions)  Let  Ti  be  a  Hilbert  space  with  inner  product  a(.,.),  such  thai 
Ji  —  'Ho  -f  •  •  •  +  Ti-n,  where  Tio, . . .  ,7in  are  subspaces.  For  v  ^  7{,  suppose  there  exists 
a  partition  v  —  Vq  -\~  •  •  •  -\-  v^,  with  v^  £  7^, ,  satisfying 


Y^a{t\,Vi)  <  n^a{v,v), 
»=o 

for  a  positive  constant  p.   Then 

n  n 

a{v,v)  <  /x^^a(P,v,P.7;)  =  p^^a(P,t;,v)  =  p'^a{Pv,v),  (2.4) 

»=0  »=0 

where  P  =  Pq  +  •  •  •  -\-  Pn,  and  the  Pi  are  the  orthogonal  projections  onto  JU  in  the 
a(.,.)  inner  product. 

PROOF  OF  LEMMA.  Let  v  =  Y^"^qI\  be  a  decomposition  satisfying  the  above  hypothe- 
sis. Then 

n  n 

a{v,v)  =  Y^a{v„v)  =  ^a(P.v,r,). 

ti=0  i  =  0 
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By  using  the  Cauchy-Schwarz  inequality,  we  obtain 

aiv,v)  <  {f;^a{P.v,P,v)y^'{f:a{vi,Vi)y^\ 

izzO  »=0 

By  using  the  hypothesis,  we  obtain 

i=0 

Squaring  both  sides  and  cancelling  the  common  terms  gives  us 

n 

a{v,v)  <  fi^J2a{PiV,v)  =  y?a{Pv,v).0 

i=0 

The  next  Lemma  introduces  sufficient  conditions  for  convergence  of  the  multiplicative 
Schwarz  method.  However,  no  estimate  for  the  convergence  factor  p  is  given. 

Lemma  16   Lei  Ti  be  a  finite  dimensional  Hilhert  space  with  inner  product  a{.,.)  and 
subspaces  ?^,  thai  sum  to  Ti: 

Then: 

(1)  there  exists  a  positive  constant  Cq,  such  that,  for  any  u  G  7{,  there  exists  a  partition 

n 

a{u,u)<cl'£a{Pu,P,u),  (2.5) 

i  =  0 

where,   as  before,   P,    denotes  the  orthogonal  projection  onto  Ti;,   in  the  a(.,.)  inner 
product. 

(2)  there  exists  a  positive  constant  p  <\  such  that  the  error  propagation  map  for  the 
multiplicative  Schwarz  method  satisfies 

\\P,tPt,---Pr'-Po^\\=P<l- 

PROOF  OF  LEMMA.     We  follow  the  proof  given  in  Lions  [22].    Consider  the  following 
map: 

{uo, . . .  ,ur,)  e  Tio  y^  •  ■  ■  X  Tin  — >  (uq  + h  u„)  G  'W. 

Let  the  norm  on  the  product  space  be  defined  by 

i=0 
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then,  clearly,  the  addition  map  is  a  continuous  map,  and  it  is  surjective  by  hypothesis. 
Thus,  by  an  application  of  the  open  mapping  theorem,  there  exists  a  continuous  right 
inverse  from  H  to  the  product  space,  i.e.,  given  u  £  7i,  there  exists  a  partition 
u  =  I3r=o'"«)  with  Ui  G  H-i,  and  satisfying 

n 

Y^a{ui,Ui)  <  cla{u,u) 

1=0 

for  some  positive  constant  cq.  By  applying  Lemma  15,  we  obtain  equation  (2.5). 

To  prove  the  second  part,  we  first  note  that  since  each  of  the  projections  are 
orthogonal,  their  norms  are  equeJ  to  one.  Thus  p  <  1.  To  obtain  a  contradiction, 
suppose  that  /?  =  1.  Then,  since  Tiis  a.  finite  dimensional  Hilbert  space,  the  maximum 
of  the  quotient  defining  the  norm  of  Pn^n-i  ' ' '  ^i'Pq'j  ^^  attained  at  some  point  u 
on  the  unit  sphere  in  Ti.  i.e.. 

Thus,  we  must  have 

\\P^u\\  =  1  =  \\u\\  =>  Pou  =  0,  and  P^u  =  u. 

By  induction,  we  obtain  that 

IIP/  •••Po^ull  =  llP^-^t^ll  =  1  =  ||u||,=^  P,u  =  0,      Vj. 

But  this  contradicts  equation  (2.5).  Thus,  our  assumption  that  p  =  1  is  wrong.  Our 
result  now  follows.  □ 

The  following  result  concerns  the  special  case  of  the  convergence  of  the  multiplicative 
Schwarz  method  in  the  case  of  two  subspaces.  In  this  case,  we  are  able  to  estimate 
the  convergence  factor  p  in  terms  of  the  parameter  Cq,  cf.  Lions  [22]. 

Lemma  17  Let  7i  be  a  Hilbert  space  with  inner  product  a(., .)  and  subspaces  Hi  and 
7^2  that  sum  to  Ti.  Assume  that  there  exists  a  a  positive  constant  cq  such  that 

Wen,       \\u\\<co{\\Pru\\'  +  \\P2u\\y^\ 

Then, 

p  =  \\P2'-Pi'\\<{i-\y^'- 
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PROOF  OF  LEMMA.     Applying  our  hypothesis  to  u  =  P^v,  we  obtain 

||P,-^r,||  <  CodiPiA-^vir  +  mP^^vWy  =  Co\\P2PM 
By  using  the  Pythagorian  theorem,  we  obtain 

\\PM'    =  \\P2P^^vr  +  mp^^vr 


Since,  V  is  arbitrary,  the  result  follows.   □ 

Many  applications  of  the  multiplicative  Schwarz  method  involve  more  than  two 
subspaces.  The  following  result,  which  is  an  extension  of  Lemma  16  due  to  Lions  [22], 
for  the  many  subspace  multiplicative  Schwarz  method,  presents  an  estimate  for  the 
convergence  factor  p  in  terms  of  other  parameters,  such  as  co,  that  can  be  estimated. 
However,  the  result  gives  a  pessimistic  and  poorer  estimate  than  Lemma  17.  Following 
that,  we  present  a  result  for  the  symmetrised  multiplicative  Schwarz  method.  See 
Widlund  [351. 

Lemma  18  Let  Ti.  be  a  Hilbert  space  with  inner  product  a(.,.).  Let  'Ho,...^'Hn  te 
subspaces  of  Ji  thai  sum  to  Ti.  Suppose  that  there  exists  a  positive  constant  Cq  such 
that  for  all  u  £  7i,  there  exists  a  partition  {u,  6  Tii}  with 

u  =  Uo  +  •  •  ■  +  Un, 

satisfying 


^a{u,,u,)  <  cla{u,u).  (2.6) 

.  =  0 

Then,  the  convergence  factor  of  the  multiplicative  Schwarz  method 

P=\\P^---Po^\\, 
satisfies 

Here  1  <  (f„,  is  given  by 
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PROOF .  We  first  prove  the  result  for  finite  dimensional  Hilbert  spaces.  Since  the  unit 
ball  in  any  finite  dimensional  Hilbert  space  is  compact,  there  exists  an  element  u  on 
the  unit  sphere  of  7i  such  that 

p  =  \\p^-.-pM- 

Since  we  have  a  product  of  orthogonal  projections, 

0</9  <  1. 


It  follows  that. 


For  i  =  0. 


\Pi^  ■••Pou\\>  p,        for  0  <  i  <  n. 


WPo'-uW  =  l=>Pou  =  0,    P^u  =  u. 
This  is  so  because,  we  can  write  u  =  PqU  +  PqU,  and 


2       \\P^---Po^A' 
p  = 


The  numerator  is  independent  of  PqU,  and  the  denominator  can  be  reduced  if  PqU  ^  0. 
Thus,  p  can  be  increased  unless  Pqu  =  0. 
Next,  We  next  show  that  for  i  >  1 

\\pur<i-p'+2{f^\\p,u\\). 

3=0 

To  do  this,  we  consider, 

||P,-tx|l^  =  1  -  WP.-'uf  -  1  -  \\P.^iP,t,---Po^)u  +  P,-'{I-Pt,---Po^M' 

=  l-(|i/^^---Po^^ir+2a(P.^(7-P.i,---Po^)u,i^^i^i,..-Po^K)+llP.^(/-i^ia---Po^)t.||^) 
<l-||P.^---Po^tz||^  +  2||P.^...Po-^u||||/^^(/-J^ii---Po')tx|| 
<l-p^  +  2||/^^(/-i^i,...Po^)u||. 

Since 

\\P.Hl-Pt.---Po^M<\\{I-Pti---Pohl 

we  consider  a  bound  for  the  right  hand  side,  which  we  write  as, 

(7  -  P.^_,  ...P,^)v  =  {I  +  P.i,  .  •  •  P,-^Po)u  -  P,ii  •  •  •  P^^u 
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=   (/  +  Pi,  ■  ■  ■  PtPo  +  •  •  •  +  Pt,P^-2  -  Pt-,>. 

Using  the  triangle  inequality,  we  obtain 

\\{l-Pt^---Po)A<\\P-i^  +  ----^Pti---PtPou\ 


«-l 


<Eli^i^ll- 

3=0 


Thus,  our  bound  for  ||Piu|l  is 


i-l 


We  have  already  shown  that 


\\Pou\\  =  0. 


(2.7) 


Denoting  r,-  =  ||P,u||,  we  obtain  that, 

f  rf     <     1  -p2  +  2(ro  +  ---^r._i)      for  i  >  1, 
I  r^     =     0  for  2  =  0. 

Using  equation  (2.6),  equation  (2.4),  and  our  definition  of  r,,  we  obtain  that 

t  =  0 

Since  each  r,  — »  0,  as  p  — >  1,  the  right  hand  side  approaches  zero  as  p  — >  1.  This 
would  lead  to  a  contradiction,  since  the  left  handside  is  1.  Thus,  p  must  be  bounded 
below  by  a  number  less  than  one,  which  depends  on  Cq.    We  now  give  a  pessimistic 


estimate  for  p,  using  square  roots  in  equation  (2.7),  and  the  fact  that   ^1  —  p"^  is 
greater  than  1  —  /?^,  for  0  <  p  <  1: 

rl     =     0 

rl     <     1-p' 


<  1-  p^  +  2>/r^ 


<3v/n^ 


r|    <    1  -  p2  +  2(v/n^  +  ^1  _  p2  +  2v/n^)    <.. 


from  which  we  obtain  by  induction  that 


,(|) 


3<(i-pT  E[1  +  '-(^+1)], 


t=i 


and  that 


p^<l- 
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To  prove  the  result  for  infinite  dimensional  spaces,  we  use  a  sequence  of  vectors  on 
the  unit  ball  such  that  the  norm  of  the  error  propagation  map  acting  on  the  sequence 
approaches  the  limit  p  as  the  sequence  approaches  its  limit.  We  then  apply  the  same 
cirguments  as  for  the  finite  dimensional  case  and  obtain  the  result.  □ 
DEFINITION.  An  iteration  of  the  symmetrised  multiplicative  Schwarz  method  con- 
sists of  2n  +  1  fractional  steps  rather  than  the  n  -f  1  fractional  steps  for  the  standard 
multiplicative  Schwarz  method.  They  involve  defect  correction  on 

no,  /ll,  .  .  .  ,  /ln-1,  /Tn,  /Tn-1,  .  .  .  ,  rll ,  /To, 

in  that  order. 

The  error  propagation  map  associated  with  the  sj/mme<nsed multiplicative  Schwarz 
method  is 

e'"^  =  ElE^e'-       where  E^  =  P^ P^_,  •  ■  ■  P,^ P,\ 

Thus,  the  error  propagation  map  is  symmetric  in  the  a(., .)  inner  product,  and  there- 
fore, the  symmetrised  multiplicative  Schwarz  method  to  solve  problem  (2.1)  can  be 
acclerated  using  the  conjugate  gradient  method. 

The  appropriate  quantity  determining  the  rate  of  convergence  of  the  symmetrised 
multiplicative  Schwarz  method  is  the  norm  of  the  error  propagation  operator,  ||U^£^„||, 
which,  by  use  of  the  previous  lemmas,  is  less  than  one.  The  quantity  determining  the 
rate  of  convergence  of  the  preconditioned  conjugate  gradient  method  is  the  quotient 
of  extreme  eigenvalues  of  7  —  E^En  ■  An  upper  bound  for  the  spectra  of  /  —  E^En  is  1, 
since  E^E^  is  a  positive  semi-definite  map,  in  the  a(., .)  inner  product.  A  lower  bound 
for  the  spectra  of  /  -  E^E^  is  1  -  p{E^Eri)  >  0,  where  p{E^En)  <  1  is  the  spec- 
tral radius  of  E^En-  Thus,  the  condition  number  of  the  preconditioned  symmetrised 
multiplicative  Schwarz  method  is  given  by  1/(1  —  p(E^En)). 

Thus,  if  we  have  a  bound  for  the  condition  number  k(/  —  E^En),  we  obtain 

^<^"^"»  = '  -  .(i-ETE„r  (^•«) 

and  this  in  turn  gives  an  estimate  for  the  spectral  radius  of  the  unsymmetrised  mul- 
tiplicative Schwarz  method: 
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The  following  Lemma,  stated  in  Widlund  [36],  gives  an  upper  bound  for  the  condi- 
tion number  of  the  symmetrised  Schwaiz  method  involving  n+1  subspaces,  in  terms  of 
the  product  of  the  condition  numbers  of  n  2  subspace  symmetrised  Schwarz  methods. 
Let  K{7io,--.,T-in)  denote  the  condition  number  of  the  symmefhsecf  multiplicative 
Schwarz  method.  Also,  we  use  K('^j,Ei=o  "^O  ^o  denote  the  condition  number  of 
the  symmetrised  multiplicative  Schwarz  method  on  "^,{=0  "^i  using  subspaces  Tij  and 
Zi=o  ^'  respectively.  Then: 

Lemma  19    The  condition  number  for  the  symmetrised  multiplicative  Schwarz  method 
satisfies 

K[7io,n\,...,nn-i,iTn)  C2  10) 

where  the  terms  on  the  right  hand  side  involve  the  condition  numbers  of  two  subspace 
symmetrised  multiplicative  Schwarz  method  using  the  subspaces  indicated. 

The  advantage  here  is  that  we  can  obtain  explicit  bounds  for  the  condition  numbers 
of  the  two  subspace  cases.  We  can  estimate  each  term  in  the  right  hand  side  of 
equation  (2.10)  by  using  equation  (2.8)  to  estimate  the  condition  number  of  the  two 
subspace  symmetrised  methods  in  terms  of  the  spectral  radius  of  the  symmetrised 
method.  The  spectral  radius  of  the  2  subspace  symmetrised  methods  can  often,  by 
using  lemma  17,  be  shown  to  be  independent  of  mesh  parameters,  etc.  Applying 
Lemma  19,  and  equation  (2.9),  we  obtain  that 

'^n  ■  ■  ■  '^1 

where  kj  =  K{T-Cj,YliZo  ^«)-  ^^  present  applications  of  this,  later  in  this  Chapter  and 
in  the  Chapter  on  iterative  refinement  methods.  We  present  a  proof  of  Lemma  19  in 
the  special  case  of  iterative  refinement  methods,  in  the  next  Chapter. 

2.1.2      The  additive  Schwarz  method. 

The  additive  version  of  the  Schwarz  alternating  method,  as  described  by  Dryja  and 
Widlund  [13],  is  based  upon  the  multiplicative  Schwarz  method,  but  modified  so 
that  the  resulting  algorithm  is  more  parallelisable.  It  is  an  iterative  method  to  solve 
equation  (2.1),  and  it  involves  the  concurrent  solution  of  subproblems  during  each 
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The  following  Lemma,  stated  in  Widlund  [35],  gives  an  upper  bound  for  the  condi- 
tion number  of  the  symmetrised  Schw&xz  method  involving  n+l  subspaces,  in  terms  of 
the  product  of  the  condition  numbers  of  n  2  subspace  symmetrised  Schwarz  methods. 
Let  k{TCo,  ■  ■  -iTi-n)  denote  the  condition  number  of  the  symmetrised  multiplicative 
Schwarz  method.  Also,  we  use  K('Wj,X^iIo  ^>)  ^^  denote  the  condition  number  of 
the  symmein'sed  multiplicative  Schwarz  method  on  J2i=o'^i  using  subspaces  Tij  and 
IIi=o  ^»  respectively.  Then: 

Lemma  19    The  condition  number  for  the  symmetrised  multiplicative  Schwarz  method 
satisfies 

kCHo, Til,..., 7in-\, Tin)  ro  in^ 

where  the  terms  on  the  right  hand  side  involve  the  condition  numbers  of  two  subspace 
symmetrised  multiplicative  Schwarz  method  using  the  subspaces  indicated. 

The  advantage  here  is  that  we  can  obtain  explicit  bounds  for  the  condition  numbers 
of  the  two  subspace  cases.  We  can  estimate  each  term  in  the  right  hand  side  of 
equation  (2.10)  by  using  equation  (2.8)  to  estimate  the  condition  number  of  the  two 
subspace  symmetrised  methods  in  terms  of  the  spectral  radius  of  the  symmetrised 
method.  The  spectral  radius  of  the  2  subspace  symmetrised  methods  can  often,  by 
using  lemma  17,  be  shown  to  be  independent  of  mesh  parameters,  etc.  Applying 
Lemma  19,  and  equation  (2.9),  we  obtain  that 

1 


piEnY   <   1  - 


5 
K„  •  •  •  Ki 


where  kj  =  fi{Tij,J2iZo  ^0-  ^^  present  applications  of  this,  later  in  this  Chapter  and 
in  the  Chapter  on  iterative  refinement  methods.  We  present  a  proof  of  Lemma  19  in 
the  special  case  of  iterative  refinement  methods,  in  the  next  Chapter. 

2.1.2      The  additive  Schwarz  method. 

The  additive  version  of  the  Schwarz  alternating  method,  as  described  by  Dryja  and 
Widlund  [13],  is  based  upon  the  multiplicative  Schwarz  method,  but  modified  so 
that  the  resulting  algorithm  is  more  paraUelisable.  It  is  an  iterative  method  to  solve 
equation  (2.1),  and  it  involves  the  concurrent  solution  of  subproblems  during  each 
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iteration,  rather  than  the  sequential  solution  of  subproblems  during  each  iteration 
as  in  the  case  of  the  multiplicative  Schwarz  method.  The  method  is  based  upon  a 
transformation  of  the  given  problem  to  a  new,  well  conditioned  linear  system  with 
the  same  solution.  Recall  that  in  the  mxiltiplicative  Schwarz  algorithm,  each  iter- 
ation involves  either  n  +  1  or  2n  -f  1  fractional  steps,  depending  on  whether  the 
unsymmetrised  or  symmetrised  version  is  used.  These  fractional  steps  are  performed 
sequentially,  although  in  some  applications,  as  we  shall  discuss  in  the  section  on  ap- 
plications to  mixed  finite  element  discretisations  of  elliptic  Neumann  problems,  many 
fractional  steps  can  be  performed  concurrently.  Thus,  the  additive  Schwarz  method 
is  more  parallelisable  within  each  iteration.  However,  the  convergence  rate  for  the 
additive  Schwarz  method  may  be  slower  than  for  the  multiplicative  Schwarz  method, 
and  the  overall  performance  may  depend  on  the  specific  problem  being  solved,  and 
the  architecture  of  the  computer  being  used,  etc. 

We  briefly  describe  the  additive  Schwarz  method.  Define  P  by 

P  =   Po  +  Pl+...  +  Pn.  (2.11) 

Note  that  since  P  is  a  sum  of  a(.,  .)-orthogonal  projections,  P  is  symmetricin  the  a(., .) 
inner  product.  Now,  given  /(.),  we  can  compute  g  =  Pu,  where  u  is  the  solution  of 
problem  (2.1),  without  knowing  u,  by  replacing  a{u,  v)  by  f{v)  in  equation  (2.2).  If  P 
is  positive  definite  in  the  a(., , )  inner  product,  then  we  could  use  an  iterative  method, 
such  as  the  conjugate  gradient  method,  to  solve  Pu  =  g  for  u.  Note  that,  since  the 
action  of  P  is  the  sum  of  various  projections,  its  action  is  highly  parallelisable.  The 
main  abstract  result  concerning  the  additive  Schwarz  method  is  given  in  the  following 
lemma.  See  Dryja  and  Widlund  [13]. 

Lemma  20  As  before,  let  Ti  he  a  Hilhert  space  with  inner  product  a(., .),  and  let  Tii 

he  suhspaces  that  sum  to  Ti: 

We  let  Pi  denote  the  a(.,.)  orthogonal  projection  onto  7i^,  and 

P  =  Po  +  P,  +  ...  +  p^. 

Then, 

\/uen,       a{Pu,u)<{n  +  l)a{u,u). 
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Suppose  that  there  exists  a  positive  constant  Co  such  that  \/u  G  "H,  there  exists  a 
partition 

U  =  Uo  +  Ui  +  •  •  •  +  Un, 

with  Ut  G  Til  satisfying 

n 

Y,a{u„Ui)  <  cla{u,u), 

1=0 

then 

Vii  G  "W,       —a{u,u)<a{Pu,u). 

PROOF  OF  LEMMA .  By  using  the  definition  of  P  as  a  sum  of  n  +  1  a(., .)  orthogonal 
projections,  we  obtain 

Ti  n  n     • 

a{Pii,u)  =  Y.a{P,u,u)  =  ^a(P,ti,P,u)  <  X^a(tx,u)  =  {n  +  l)a{u,u). 

t  =  0  j  =  0  »=0 

We  use  our  hypothesis  and  Lemma  16  to  obtain  the  lower  bound  for  P  in  the  a(.,.) 
inner  product: 

—a{u,u)  <  a{Pu,u). 

This  shows  that  P  is  positive  definite  in  the  a(., .)  inner  product.  The  rate  of  conver- 
gence of  the  conjugate  gradient  method  to  solve  Pu  =  g,  is  bounded  in  terms  of  the 
condition  number  of  P  which  is  bounded  by  (n  -f  Ijcg.  Q 

Dryja  and  Widlund,  see  [13],  have  studied  this  algorithm  in  the  case  of  standard 
finite  element  discretizations  of  elliptic  problems,  in  2  and  3  dimensions.  They  show 
that  in  the  case  where  one  of  the  subspaces  is  a  coarse  model,  and  the  other  subspaces 
are  suitably  chosen,  then  there  exists  positive  constants  61,62,  independent  of  the 
mesh  parameters,  such  that 

61  a{u,u)  <  a{Pu,u)  <  62  a[u,u),     Vu  G  S. 

Thus,  if  the  conjugate  gradient  method  is  used  to  solve  Pu  =  g,  in  the  a(.,.)  inner 
product,  the  number  of  iterations  would  depend  on  62/61,  which  is  independent  of 
the  mesh  parameters.  Numerical  tests  run  by  Greenbaum,  Li  and  Chao  [18]  verify 
the  theoretical  results  obtained. 

We  also  mention  that  X.  C.  Cai  has  studied  a  generzJisation  of  the  additive 
Schwarz  method  to  ceitain  non-symmetric  problems,  including  convection  diffusion 
problems,  and  time-dependent  diffusion  problems.  See  Cai  [11]. 
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2.1.3      The  Schwarz  methods  for  saddle  point  problems. 

We  now  consider  the  steps  involved  in  applying  iterative  methods  to  saddle  point 
problems.  In  particular,  we  consider  the  multiplicative  and  additive  Schwarz  methods. 
See  also  Lions  [22],  and  Glowinski  and  Wheeler  [16].    Consider  the  discrete  saddle 
point  problem  obtained  by  discretising  the  continuous  problem  (1.7). 
Find  Uh  E  Xh,Ph  ^  Yh  such  that: 

(Ah    B^\  (  uh\  _  (  Wh 


J^h      ^     J   \  Ph  J        \   rh   J 

where  Ah  is  a  symmetric,  positive  definite  matrix,  and  Bh  is  of  full  rank.  With  some 
abuse  of  notation,  we  use  Uh,Ph,  etc.,  to  denote  functions  as  well  as  their  vector 
representation  in  matrix  notation. 

Note  that,  for  discretisation  of  elliptic  Neumann  boundary  value  problems,  ph  is 
unique  only  in  A''/,(Q)  =  Qh{^)/R,  i.e.,  unique  up  to  a  constant  in  Qh{Q).  Thus,  if 
we  use  a  basis  for  Qh{^)  in  our  construction  of  the  matrix,  then  B^  will  have  a  null 
space  spanned  by  (1, . . . ,  1)-^.  We  refer  to  functions  in  Xh  as  velocities,  and  functions 
in  Yh  as  pressures;  B^  will  be  refered  to  as  the  discrete  divergence  map. 

Let 

n  =  {uHe  Xh  :  BhUh  =  0}. 

Elements  in  7i  will  be  denoted  with  a  subscript  DF,  for  example  Vh,DF  is  discrete 
divergence  free.  We  can  decompose,  in  many  ways,  the  velocity  space  ■ 

Xh  =  Xh,i  +  Ti; 

Xh,i  can  even  be  taken  to  be  A^.  Let  A'o, . . . ,  A'n  be  subspaces  of  Xh  that  sum  to  Xh- 
Let  Yo,...,Yn  be  the  ranges  of  the  Bh  map  acting  on  A''o,...,A'n  respectively.  Let, 
Tii  denote 

n,  =  Xinn  =  {vh  €  Xi  :  qJ^BhVh  =  0,     Vg;.  G  Vj}. 

Thus,  Tii  also  sum  to  7i.  We  divide  the  solution  process  into  three  steps. 

1.  In  the  first  step  we  compute  a  velocity  Uhj  G  Xh,i  such  that  BhUhj  =  Fh  =  BhUh- 
In  applications,  we  usually  solve  a  small  coarse  mesh  problem  followed  by  the  con- 
current solutions  of  n  subproblems  to  determine  u/,  /. 
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2.  In  the  second  step,  we  compute  the  discrete  divergence  free  velocity 

Uh,DF  =  Uh  -  Uh.i  €  v.. 

Note  that,  {uh,DF,Ph)  satisfies  the  following  equation: 

(Ah    J5M  /  Uh,DF  \  ^(  Wh-  Anuhj  \  ^(  ^H-  AhUHj  \ 
[  Bh      0        [     PH  [  Fh-  Bhuhj  [  0 


(2.13) 


Thus,  by  restricting  the  problem  to  discrete  divergence  free  velocities,  we  find  that  if 
(uh,DF,Ph),{vh,DF,qh)  eH  xYh,  then 

/  Vh,DF   \       (Ah      B^   \    (   Uh,DF   \   _    (  ■Vh,DF  \       (   Wh-  AhUhJ  \  (o^A\ 

[      qh      )      [Bh       0     )[      PH      )-[      q,      )      [  0  j  ■  ^^-^^^ 

This  leads  to  an  equivalent  problem  to  determine  Uh,DF'- 

{Find  Uh,DF  £  "^j  such  that  .        . 

Vh.DF'^hUh.DF  =  l^h.DFi^^^h  -  AhUhj),         Vv/,  G  7^. 

If  Ah  is  positive  definite  on  Ti,  it  follows  that  the  problem  to  find  Uh,DF  is  symmetric, 
positive  definite.  We  use  iterative  methods,  such  as  the  Schwarz  methods  based  on 
the  subspaces  7i,  to  solve  for  Uh,DF  in  equation  (2.15)  Note  that,  standard  iterative 
methods  to  solve  equation  (2.15),  may  require  that  the  problem  be  expressed  in  terms 
of  a  basis  for  Ti  and  that,  such  a  basis  may  be  computationally  expensive  to  construct. 
However,  we  can  implement  the  multiplicative  and  additive  versions  of  the  Schwarz 
alternating  method  without  the  use  of  a  basis  for  the  discrete  divergence  free  velocity 
space  Ti.  We  find  the  projections  onto  Tii  by  solving  a  saddle  point  subproblem  on 
Xi  X  Yi-  The  second  step  is  the  most  expensive  step  in  the  solution  process.  In  most 
applications  involving  the  multiplicative  and  additive  Schwarz  methods,  the  number 
of  iterations  required  are  independent  of  the  mesh  width  h.  During  each  iteration, 
many  subproblems  can  often  be  solved  in  parallel. 

3.  Finally,  we  compute  the  pressure  ph-  In  most  applications,  the  computational 
cost  of  this  step  is  the  same  as  or  less  than  the  cost  of  one  iteration  in  step  2.  One 
could  use  various  subproblems  to  determine  pressures  in  subspaces,  and  then  compute 
the  global  pressure  ph-  We  discuss  the  specific  details  later  on  in  this  chapter  and  in 
the  chapter  on  iterative  refinement  methods. 
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In  the  next  section,  we  apply  these  methods  to  the  particular  case  of  mixed  fi- 
nite element  discretisations  of  elliptic  Neumann  problems  using  the  Raviart-Thomas 
finite  element  spaces.  Recall  that  the  Raviart-Thomas  finite  element  spaces  have  the 
property  that 

discrete  divergence  free  <=^  L^- divergence  free. 

2.2  Application  of  the  Schwarz  methods  to  mixed 
finite  element  discretisations  of  elliptic  Neu- 
mann problems  using  the  Raviart-Thomas  fi- 
nite element  spaces. 

Let  f]  be  a  polygonal  region  in  R^  triangulated  by  a  uniformly  shape  regular  trian- 
gulation  r^  which  is  a  refinement  of  a  coarser  triangulation  r^  of  width  H.  Thus, 
every  element  in  t"  can  be  expressed  as  a  union  of  elements  in  t^.  The  triangula- 
tion consists  entirely  of  triangles  or  entirely  of  parallelograms.  Let  the  elements  in 
the  coarse  triangulation  r^,  or  as  we  will  refer  to  them,  subdomains,  be  denoted  by 
Qi,...,QAr.  For  each  subdomain  S7,,  we  construct  an  extension  Q"'  with  f2,  C  Q"*. 
The  extension  fi"'  is  assumed  to  have  (ii6ia7ice(5r2,5f2"')  >  cH,  for  some  positive 
constant  c,  independent  of  i  and  H.  In  case  f]"'  C  Cl,  we  define 

In  case  Q"'  ^  Q,  we  cut  off  the  portion  of  J7"'  not  in  J7,  i.e., 

See  Figure'2.2.  We  assume  that  the  extended  subdomains  0,'-  are  unions  of  elements 
in  the  triangulation  r'*  of  Q. 

For  a  fixed  non-negative  integer  r,  let  {Xh{^),Yh{^))  denote  the  coarse  mesh 
Raviart-Thomas  space  of  order  r  for  the  Neumann  problem;  i.e., 

and 

YH{n)  =  QHin)r\L\n)/R. 
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Similarly,  {Xh{^),  Yhi^))  denotes  the  Raviart-Thomas  spaces  of  order  r,  based  on  the 
fine  triangulation  r'*,  for  the  Neumann  problem.  Let  Ijj  :  Xh{^)  — >  Xh{^)  denote 
the  interpolation  operator  for  the  velocity  space.  Since  Xh{^)  C  Xh{^),  Ig  is  the 
identity  on  Xh{^)-  As  before,  with  some  abuse  of  notation,  we  use  Uh,UH,  Ph,PH,  etc., 
to  denote  functions  in  the  Raviart-Thomas  finite  element  spaces,  and  interchangeably 
their  representation  with  respect  to  the  basis  used  in  in  the  matrix  notation.  Their 
use  will  be  clear  from  the  context. 

As  discussed  in  the  previous  section,  the  solution  of  the  saddle  point  problem 
consists  of  three  steps.  Step  1  and  step  3  are  the  same  for  both  the  miiltiplicative 
and  additive  Schwarz  methods.  In  the  case  of  the  multiplicative  Schwarz  method, 
the  pressure  ph  can  be  determined  in  step  3  without  the  use  of  subproblem  solves, 
and  may  involve  only  on  the  order  of  N  floating  point  operations.  For  step  2,  we 
prove  that  the  rates  of  convergence  of  both  the  multiplicative  and  additive  Schwarz 
methods  are  independent  of  h.  Numerical  results  are  presented  in  a  later  section. 
They  indicate  that  the  rates  of  convergence  for  both  the  multiplicative  and  additive 
Schwarz  methods  are  independent  of  h  and  H,  the  fine  and  coarse  mesh  parameters, 
respectively,  but  depend  mildly  on  c  the  amount  of  overlap  between  subdomains. 

2.2.1      Step  1:  Reducing  the  saddle  point  problem  to  a  sym- 
metric positive  definite  problem. 

To  reduce  the  saddle  point  problem  to  a  symmetric  positive  definite  problem,  we  need 
to  compute  a  discrete  velocity  Uh,i  €  Xh{^)  satisfying: 

BhUk.i  =  Fh  -  BhUh- 

First,  we  solve  a  coarse-mesh  problem  to  obtain  a  coarse  mesh  discrete  velocity 

uh  with  normal  trace  on  dVtj  compatible  with  the  mean  value  of  /(x)  on  Q,j  for  each 

;,  i.e., 

Find  UH  G  Xh{^),  Ph  e  Yh{^)  such  that 
a{uH,VH)  +  b{vH,PH)     =    0  ^Vff  e  Xn{^)  (2.16) 

b{uH,qH)  =    F{qH)  =  Inf{x)qndx     "^qn  e  Ynin) 

Substituting  qn  =  Xn,  €  Yh{0,),  the  characteristic  function  of  fij,  into  the  divergence 

constraint  satisfied  by  uh,  we  see  that: 

/    f{x)dx  =   f   V-undx^  I     n-unds,  (2.17) 

Jn,  Jn.  Jan. 
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i.e.,  ujj  has  a  normal  trace  on  dCli  which  is  compatible  with  the  mean  value  of  f{x) 
on  Cli  ioi  i  =  I,. ..  ,N.  Now,  u^  —  I^un  has  the  property  that 

has  mean  value  zero  in  each  of  the  subdomains  fii, . . . ,  Vlj^.  Thus,  if  we  assign  zero 
normal  trace  on  each  subdomain  boundary  5^2,,  we  would  have  compatible  data  to 
pose  N  parallel  subproblems  on  fii , . . . ,  Q,i^  for  local  velocities  and  pressures  {6uf^,Sp\) 
using  the  new  right  hand  side,  given  by  the  residual,  as  follows: 

Find  5u\  e  Xhi^i),  Sp[  G  y^(J7.)  such  that 
a{Sul,v)  +  b{v,8p'^)     =     -a{I^un,v),  W  e  Xh{^i) 

[     HSllq)  =     Fiq)-b{I^UH,q),      Vg  G  ^(a). 

The  local  spaces  (A'"h(r2,),  ^/^(n,))  are  defined  by 

Y,{n,)^QH{n)nL'{no/R- 

Since  the  local  problems  have  compatible  boundary  conditions,  see  (2.17),  they  are 
well  posed.  Since  the  normal  traces  are  zero  on  interface  boundaries,  they  match;  thus 
this  composite  velocity  will  be  in  A';,(fi).   See  Lemma  5  about  composite  functions. 
We  form  the  composite  velocity  by  adding  up  the  subdomain  velocities  6uf^. 
Let 

By  construction,  Uh,!  satisfies  BhUhj  —  BhUh  =  F^.  Thus  Uh,DF  =  ^h  —  'Uh,i,  satisfies: 

j  a{uh,DF,v)  +  b{v,Ph)    =     -a{uh,i,v),      W  e  Xh{^)  (o -iR) 

1  b{uH.DF,q)  =     0,  Vgen(fi),  ^      ^^ 

i.e.,  Uh,DF  is  divergence  free.  We  now  show  that  finding  Uh,DF  is  a  positive  definite 
symmetric  problem  on  the  divergence  free  subspace  7i  defined  by: 

n  =  Vf,{n)  n  Ho,an{div\  fi).  (2.19) 

From  equation  (2.18),  we  see  that  Uh,DF  £  "H.  Restricting  v  in  problem  (2.18)  to  Ti, 
we  obtctin  the  following  equivalent  problem: 

Find  Uh.DF  €  'H  such  that 

a{-Uh.DF,v)  -  -a{uh,i,v)        Vv  G  7i.  ^'     ' 
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Since 

and  since  a(.,.)  is  ^((fu'°,Q)-elliptic  by  equation  (1.13),  it  follows  that 

When,  a{^h,Vh)    >    Cx\\Vh\\)}^j-^f^y 

Therefore,  problem  (2.20)  is  a  symmetric  positive  definite  problem. 

2.2.2      Step  2:     Solution   of  the  divergence  free  symmetric 
positive  definite  problem. 

We  solve  problem  (2.20)  by  the  use  of  the  multiplicative  or  additive  Schwarz  methods. 
Once  Uh,DF  is  obtained,  the  discrete  velocity  solution  Uh  is  determined  by  letting: 

Uh  =  UhJ  +  Uh,DF- 

First,  we  introduce  some  notation.  For  i  =  1, . . . ,  iV,  let: 

X^m    =    Vk(Q)n^o,an'(^it',f^:)  /,  21) 

These  are  subspaces  of  the  Raviart-Thomas  spaces  obtained  by  restricting  them  to 
the  extended  subdomains  {Cl[}.  The  functions  in  these  subspaces,  when  extended  by 
zero  outside  Q'-  to  the  whole  domain  Cl,  lies  in  the  global  Raviart-Thomas  spaces. 

The  spaces  used  in  the  formulation  of  the  multiplicative  and  additive  Schwarz 
methods  are  the  divergence  free  subspaces  of  the  Raviart-Thomas  velocity  spaces, 
i.e., 

-Ho  =  VH{n)nHoM^iv°,n)  (2.22) 

n,  =  VH{n)nHo,au'Sdiv°,Q[). 
As  stated  before,  we  do  not  possess  a  basis  for  the  divergence  free  subspaces 
7i,7io.,. . .  jTiN-  To  implement  the  Schwarz  methods,  we  need  to  compute  the  pro- 
jections onto  the  divergence  free  subspaces.  For  i  =  0,...,iV,  let  Pi  denote  the 
orthogonal  projection  in  the  a(.,.)  inner  product  onto  Hi.  We  sometimes  use  Ph  to 
denote  Pq,  the  orthogonal  projection  onto  the  coarse  mesh  divergence  free  velocity 
space  Hq.   To  obtain  each  projection  we  solve  a  saddle  point  subproblem  using  the 
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basis  for  the  appropriate  Raviart-Thomas  subspaces.  For  instance,  to  find  PiW,  for 
i  =  1, . . . ,  A'^,  we  solve: 

Find  WH  e  XH{n[),  SH  e  n(r2;)  s.t 

a{wH,v)  +  b{v,SH)  =     a{w,v)      Wv  e  Xhi^^)  (2.23) 

b{wH.q)  .  =     0  VgelUn;.) 

Thus  P,w  =  whe  Tii. 

For  N  sufficiently  large,  let  us  partition  the  collection  of  subdomains  fii,...,fi';v 

into  n  colors  1, . . . ,  n,  coloring  two  subdomains  the  sanae  only  if  they  are  disjoint.  Let 

J,-  denote  the  collection  of  indices  of  the  subdomains  {Q'-}  of  color  i,  and  let  d  be 

the  union  of  all  subdomains  {0.'^}  of  color  i.  Thus, 

C,  =  U,•eJ.^;.. 
Fo^  all  shape  regular  triangulations  t^ ,  the  number  of  colors  n  is  independent  of  H 
and  h.  One  can  show  this  using  the  fact  that  the  number  of  neighbouring  elements 
for  a  shape  regular  triangulation  is  bounded  independently  of  the  mesh  widths.  For 
instance,  in  the  case  of  a  rectangular  region  triangulated  by  an  uniform  mesh  of 
rectangles,  the  number  of  colors,  for  suitably  chosen  overlap  among  the  extended 
subdomains,  is  four,  independent  of  H  and  h.  See  Figure  2.3,  in  the  section  on 
Numerical  results. 
Define 

T^c.  =  i:^.'  (2.24) 

for  each  color  i  =  I,...,??.  This  is  an  orthogonal  direct  sum.  We  now  state  and  prove  a 
result  concerning  the  existence  of  a  partition  for  the  divergence  free  Raviart-Thomas 
subspaces.  This  result  will  be  used  to  prove  the  basic  convergence  results  for  the 
multiplicative  and  additive  Schwarz  methods  applied  to  elliptic  problems  discretised 
by  the  Raviart-Thomas  mixed  finite  element  method. 

Lemma  21  Let  T^c'Wci  •■•7^c„  ^e  subspaces  of  TC  as  defined  by  equation  (2.24)- 
Then,  there  exists  a  positive  constant  /i,  independent  of  h  such  that  for  every  Vh  G  Ti, 
there  exists  a  partition  with  Hi  G  "Wc, ,  for  i  =  1, .  . .  ,n,  and  Vq  G  "Ho  satisfying 

Vh  =  Vo  +  Vi  + h  v„, 

and 

n 

^a(f,,v,)  <  fi^a{ijh,Vh). 

i  =  0 
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Figure  2.1:  Partitioning  Ci  into  0,  and  A,-. 


PROOF  OF  LEMMA.     First,  we  define  vq  =  Pqv.  Since  Pq  is  a  projection,  we  obtain: 

(2.25) 


a{v  —  Vo,i^  —  vo)  <  a{v,v). 


Next,  we  construct  a  partition  for  v  —  vq  =  Vi  +  . . .  +  r„,  with  each  r,-  G  Tic,  having 
support  in  C,.  We  construct  Vi,...,Vn  sequentially,  starting  with  Vi.  To  describe 
the  construction  of  {T,  at  step  z,  we  partition  the  colored  region  C,-  into  two  disjoint 
sets  0,-  and  A,.  Though  the  region  C,  may  not  be  connected,  it  is  a  union  of  one  or 
more  of  the  connected  subdomains  {^'j}.  The  subregions  0,-  and  A^  will  also  be  the 
union  of  connected  sets.  0;  is  defined  as  the  subset  of  C,-  which  does  not  intersect  the 
remaining  regions  C,+i, . . . ,  C„.  A,  is  defined  as  the  complement  of  0^  in  C,-.  Thus, 
for  i  =  1, . .  .  ,n: 

\  A.     =    C.  -  0.  =  Ci  n  0,^ 
See  Figure  2.1.  Note  that  0„  =  C„,  A„  =  0,  and  U^^i©.  =  H. 

Once  Vo, . . . ,  Vi_i  have  been  defined  on  UjljCj,  we  define  Vi  on  0j.    Since  the 
remaining  regions  A,,  Ci+i, . . . ,  Cn  do  not  intersect  0^,  we  have  no  alternative  but  to 

define 

t-i 

(2.26) 


t7,  =  v  —  2J  Vj  on  0, . 

j=0 
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Since  the  Vj  are  defined  to  be  zero  outside  C^,  for  jf  =  z  +  1,  •  •  •  ,n,  this  gives  us 

Vo  +  vi  +  .. .  +  Vn  =  V      on  0.. 

By  the  definition  of  v,  on  0,  the  norm  of  Vi  on  0,  equals  the  norm  of  i7  —  Yi)Z])  ^j  O" 
Qi.  Thus, 

Next,  we  describe  how  {T,  is  defined  on  A,.  Since  Vi  is  to  have  support  in  Ci,  the 
normal  trace  of  Vi  on  5C,  must  be  zero.  We  must  also  make  the  normal  trace  of  v,-  on 
(9A,  compatible  with  the  normal  trace  of  {T,  on  50,,  wherever  5A,  and  dQi  intersect; 
Vi  is  defined  on  0,  by  equation  (2.26).  Therefore,  if  we  let  ^,  denote  the  normal  trace 
n  •  Vi  on  5A,,  we  must  have: 

/  5.  =  -ln[v  -  L}=o  ^i\      on  5A,  n  50,,  ,         . 

\  5,  =  0  on  5A,  n  5C.  =  5A,  -  (50.  n  5Ai).  ^       ' 

We  claim  that  Jq,^  Qidsx  =  0.  There  are  two  cases  to  consider.  To  simplify  our 
proof,  we  consider  each  connected  subdomain  Q'-  that  constitutes  C,-  separately.  See 
Figure  2.1. 

1.  In  case  0.^^  C  fi,  then  fi^  n  A,  will  be  ring  shaped,  and  dQi  fl  0.'-  will  be 
5Ai  n  dQi  n  r2' .  Thus,  since  gi  is  the  normal  trace  on  dQi  fl  fi'  of  i7  —  X^'Cq  Vj, 
which  is  divergence  free,  our  claim  follows. 

2.  In  case  f2"'  ^  Q,  then  Q^nA,-  will  be  'L'  or  'U'  shaped  with  part  of  its  boundary 
intersecting  dVt.  Since  v  —  Yl^l^o'^j  ^^^  2^^°  normal  trace  on  5f2  and  since 
■i;  —  Z!j=o  ^j  ^s  divergence  free  on  0,  fl  Jl^-,  it  follows  that  /^^  nae  nfi'  Qi'^^x  =  0. 

From  this,  our  claim  follows.  Furthermore,  by  an  application  of  the  normal  trace 
lemma,  we  can  bound 

II^.IIh-./.oa.)  <  c(A.)||i;.i|^(,,„e.)-  (2-28) 

Using  the  normal  trace  g,  of  v^  defined  by  equation  (2.27)  on  5Ai,  we  define 

ViU,  =  E^gi       on  A„  (2.29) 

where  E  denotes  the  Raviart-Thomas  divergence  free  velocity  extension  of  the  nor- 
mal trace  g,.  Such  an  extension  exists  by  the  extension  theorem  for  Raviart-Thomas 
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velocity  spaces  described  in  Chapter  1.  Furthermore,  the  extension  theorem  guaran- 
tees that  there  exists  a  positive  constant  c  =  c(Ai)  with 

ll^.|lH(ci..,A.)  =  ll^'5.b(d.„,A.)   <  <^(A0ll5.||i7->/:OA.)-  (2-30) 

Combining  equations  (2.26),  (2.28),  (2.29),  (2.30),  we  obtain,  for  some  positive 
constant  c(Ai,  0,) 

II^"'.|Ih(..v,c.)  ^  ^(A.'0.)(E  IWM'i'v.o))-  (2.31) 

i=o 

By  using  equation  (2.25),  equation  (2.31),  and  induction,  it  follows  that  there  exists  a 
positive  contant  c,,  depending  on  the  various  regions,  and  n,  which  are  all  independent 
of  h,  such  that: 

a{vt,Vi)  <  c^a{v,v). 

Finally,  summing  over  2,  we  obtain,  for  some  positive  fi  independent  of  h,  the  follow- 
ing: 

n 

REMARK.  The  result  holds  without  the  use  of  ?^o  and  Vq  in  the  partition,  i.e.,  given 
V  ^  7i,  there  exists  a  partition 

V   =  V^   + \-Vri, 

with  Vi  G  Ti-c,-,  satisfying 

n 

^a(i7,,{ri)  < /i^a(r,{r), 

t=i 

for  some  positive  constant  fi  independent  of  h  and  v  ^  7i.  The  proof  is  the  same. 

Using  Lemma  21,  we  prove  the  main  results  about  the  rate  of  convergence  of 
the  multiplicative  and  additive  Schwarz  methods  applied  to  the  mixed  finite  element 
discretisation  of  elliptic  problems.  First,  we  prove  the  following  result  for  the  unsym- 
metrised  multiplicative  Schwarz  method. 

Theorem  3   For  the  multiplicative  Schwarz  method  using  the  subspaces  7io,7ici,  •••  ■,'Hc„ 
or  the  subspaces  Tici  ■>■■■■>  ^c„;  ^he  convergence  factor  p  is  bounded  less  than  one,  in- 
dependent of  h. 
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PROOF.     The  result  follows  by  use  of  Lemma  21  and  Lemma  18,  since  the  constants 

involved  are  inedpendent  of  /z.  □ 

The  following  result  concerns  the  additive  Schwarz  method. 

Theorem  4  P  =  Pq  +  Pi  +  •  •  •  +  Pn  is  symmetric  and  positive  definite  in  the  a(., .) 
inner  product.  Furthermore,  there  exists  positive  constants  61,62,  independent  of  h, 
such  that 

61  a{uh,Uh)  <  a{Puh,Uh)  <  61  a{uh,Uh)- 

PROOF  OF  UPPER  BOUND.  P  is  symmetric  since  it  is  a  sum  of  orthogonal  projections, 
each  of  which  is  symmetric,  in  the  a(.,.)  inner  product.  Using  the  coloring  of  the 
subdomains  as  before,  Ci, . . . ,  Cn,  we  may  define  for  each  z, 

Pc.  ^  E  P.- 

Each  Pc,  is  an  orthogonal  projection,  since  it  is  a  sum  of  orthogonal  projections  onto 

subspaces  that  are  mutually  orthogonal.    Then  P  =  Pq  +  Pc^  +  ...  +  Pc„,  and  so 

the  upper  bound  ^''2  is  n  +  1.    Note  that,  n,  the  number  of  colors  for  shape  regular 

triangulations,  is  bounded  independently  of  H  and  h.  □ 

PROOF  OF   LOWER  BOUND  .     For  the  lower  bound,  we  use  Lemma  15  and  Lemma  21  to 

obtain: 

-^a{v,v)  <  a{Pv,v), 

with  a  constant  fi  independent  of  h,  but  possibly  dependent  on  the  geometry  of  the 
subdomains. 

2.2.3      Step  3:   Determining  the  pressure. 

We  determine  the  pressure  p/,  as  follows.  Let  (uh,ph)  denote  the  exact  solution  of  the 
discrete  problem  (1.14).  Then,  by  construction, 

Uh  =  Uh,I  +  Uh,DF, 

which  was  computed  in  the  previous  two  steps.  Recall  that,  for  the  mixed  formulation 
of  elliptic  Neumann  problems,  Wh  =  0,  thus,  ph  satisfies  the  following  equation: 

Find  Ph  G  Yh{n)  s.t 

a{uh,i  +  Uh,DF,v)  +  b{v,Ph)     =    0      WeXh{n)  (2.32) 

b{uh.i  +  Uh,DF,q)  =     0      \/qeYhin). 
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Now,  we  can  find  the  projection  onto  the  divergence  free  Raviart-Thomas  subspaces, 
for  each  z,  using  the  new  residual,  i.e.,  for  z  =  1, . . . ,  A'^  solve: 

'  Find  w,^H  e  Xh{Cl[),p,_h  £  Yhi^d,  such  that 

*      a{w,^h,v)  +  b{v,pi^h)     =    0  -  a{uhj +  Uh,DF,v)  =  K'^,Ph),     VtT  G  J\:;,(Q-) 
,     Kw,_f,,q)  =    0,  Vgen(n:.) 

(2.33) 

Each  Wi^h  =  0,  but  pih  will  be  nonzero  in  general.  Since  the  divergence  operator  maps 

X/,(Q()  onto  Yhi^'i),  'w^  obtain  that 


/   {p,_H  -  Ph)qdx  =  0,       \/qeY,,{n[). 


i.e,  pi^h  is  the  L^  projection  of  p/j  onto  y/,(QJ).  By  using  the  uniqueness  up  to  constants 
of  the  pressure  solutions,  we  obtain,  for  2  =  1,...,  N ,  that  : 

Pi,h  =  Ph  +  c,     in  f}-, 

for  some  constants  c,.  Thus,  in  the  intersection  between  the  subdomains,  we  must 
have 

p,.h  =  P:.h  +  Ci  -  Cj      in  Cl'-  H  Q'j. 

Thus,  once  {p,,h}  are  computed,  we  can  add  an  appropriate  constant  in  each  extended 
subdomain  fi(  for  each  i  =  1,...,A'^,  sequentially  to  get  a  consistent  pressure  in 
U'-jO'.  In  the  last  step  we  obtain  a  global  pressure  which  differs  from  the  individual 
subdomain  pressures  {pi^h}  by  some  constant  in  each  subdomain  Q,[.  When  this  global 
pressure  is  normalised  to  have  mean  value  zero  in  Q,  it  equals  the  exact  pressure 
solution  ph  defined  by  equation  (2.12). 

2.3      Numerical  results. 

In  this  section,  we  present  the  results  of  numerical  tests  with  the  multiplicative  and 
additive  Schwarz  algorithms  for  mixed  finite  element  discretisations  of  elliptic  prob- 
lems. We  choose  the  domain  of  the  elliptic  problem  to  be  the  unit  square  [0, 1]  x  [0, 1] 
in  the  x  —  y  plane  and  let  the  coefficients  A{x,y)  of  the  elliptic  operator  be  diagonal 
and  piecewise  constant: 

A{x,y)  =  a{x,y)I, 
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where  /  denotes  the  identity  matrix  and  a{x,y)  is  a  piecewise  constant  function  with 

a  jump  across  x  =  |  : 

/  1     0  <  X  <  0.5 
"^^'^^=1  J    0.5<x<l  (2-^^^ 

The  ratio  of  the  coefficients,  J,  varies  between  1  and  10~®.  To  obtain  the  discretisa- 
tion, we  use  the  lowest  order  Raviart-Thomas  spaces  defined  on  a  rectangular  mesh. 
We  note  that  for  a  uniform  mesh,  such  as  in  Figure  2.3,  we  can  represent  the  Raviart- 
Thomas  spaces,  as  a  tensor  product  of  two  spaces,  one  along  each  axis,  of.  Glowinski 
and  Wheeler  [16]. 

For  a  positive  integer  n,  the  fine  mesh  on  the  unit  square  is  taken  to  be  the  n  *  n 
uniform  square  mesh.  Let  h  =  1/n  denote  the  fine  mesh  parameter,  then  the  diameter 
of  the  square  element  is  %/2/i.  The  coarse  mesh  is  obtained  by  dividing  the  square  into 
n,  *n,  square  elements,  called  subdomains.  We  assume  that  n,  divides  n.  Thus,  each 
subdomain  (coarse  mesh  element)  Q,  contains  (n/n,)*  (n/n,)  fine  mesh  elements.  We 
let  the  coarse  mesh  parameter  he  H  =  l/n,;  thus  the  diameter  of  the  coarse  mesh 
element  is  \/2/n,.  To  obtain  an  overlap  between  the  subdomains,  we  introduce  an 
integer  parameter  Hq.  The  extended  subdomains  are  obtained  by  centering  a  square 
Q"'  of  (TJo  +  n/rij  -f  TZo)  *  {jio-^n/n,  +  no)  elements  at  the  center  of  subdomain  fi,-,  and 
removing  the  part  outside  the  unit  square.  The  resulting  subregions  are  denoted  by 
QJ  =  Q"'  n  Q.  See  Figure  2.2,  for  an  example  without  the  elements  in  the  fine  mesh 
drawn  in.  Thus,  in  the  case  of  the  interior  subdomains,  Q,[  contains  {n/n,  +  272,,)^ 
elements.  The  results  are  presented  for  various  choices  of  n,n,,no  and  for  various 
choices  of  discontinuities  in  the  coefficients  of  the  elliptic  operator. 

We  make  a  few  remarks  about  the  linear  system  obtained  by  using  the  lowest  order 
Raviart-Thomas  finite  element  spaces  on  a  uniform  rectangular  mesh.  We  recall  that 
the  linear  system  has  the  form 


Bh      0 


Vh 


Fh 


(2.35) 


where  Ah  is  a  symmetric,  positive  definite  matrix.  Since  the  velocity  unknowns  are 
the  values  of  the  normal  component  on  the  edges  of  the  rectangular  elements,  one 
on  each  edge,  we  have  as  many  velocity  unknowns  as  there  are  element  edges  in  the 
interior  of  the  unit  square.  Thus,  there  are  2n(n  —  1)  such  degrees  of  freedom.  For 
diagonal  A{x,y),  the  velocity  unknowns  on  the  vertical  edges  are  uncoupled  from  the 
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Figure  2.2:  Su 


adomains  and  extended  subdomains. 
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velocity  unknowns  on  the  horizontal  edges  in  the  matrix  Ah]  however,  in  the  stiffness 
matrix,  there  will  be  couplings  via  the  pressure.  For  a  suitable  ordering  of  the  velocity 
unknowns,  the  leading  block  Ah  of  the  stiffness  matrix  is  tridiagonal,  with  couplings 
between  adjacent  vertical  edges  and  adjacent  horizontal  edges. 

For  the  Laplacian,  i.e.,  A{x,y)  —  I,  one  can  verify  that  Ah  is  well  conditioned, 
uniformly  in  h.  However,  if  A{x,y)  contains  a  discontinuity,  then  the  condition 
number  of  Ah  will  be  proportional  to  the  magnitude  of  the  jump,  uniformly  in  /i,  and 
for  large  discontinuities,  the  matrix  Ah  can  be  ill  conditioned. 

If  a  basis  for  the  pressure  space  Qh{^)  is  used  in  the  construction  of  the  stiffness 
matrix,  instead  of  Ih(^)  =  Qh{^)/R^  then  the  stiffness  matrix  has  a  null  space 
corresponding  to  a  constant  pressure.  Thus,  in  that  case,  Bh  has  a  null  space  of 
dimension  1  corresponding  to  a  constant  pressure.  Recall  that  for  the  lowest  order 
Raviart-Thomas  space  on  a  rectangular  mesh,  there  is  one  pressure  node  associated 
with  each  rectangular  element.  Thus,  the  dimension  of  Qh{^)  is  n^-  We  note  that 
the  block  matrix  Bh  does  not  contain  any  information  about  the  coefficients  of  the 
elliptic  operator,  A{x,y),  since  it  is  an  approximation  of  the  divergence  map. 

The  Schur  complement  associated  with  ph  is  BhAh^B^,  and  it  is  ill  conditioned 
as  h  — >  0,  even  for  A{x,y)  =  I.  This  is  because  this  Schur  complement  represents  a 
discretisation  of  the  elliptic  problem  in  p.  We  can  show  that  the  Schur  complement 
is  ill-conditioned  by  using  the  tnf-sup  condition  and  the  fact  that  the  a(.,.)  norm  is 
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not  equivalent  to  the  H{divM)  norm,  except  for  divergence  free  functions. 

Next,  we  discuss  a  preconditioned  iterative  method  that  we  use  to  solve  the  local 
problems.  Let 

(2.36) 


'Ah    Br 
.Bh      0 

= 

denote  a  locaJ  stiffness  matrix.  Note  that,  the  extended  subdomains  in  the  interior 
contain  2(n/n,  +  2no){n/n,  +  2no  —  1)  velocity  unknowns,  and  {nin,  +  2noY  pressure 
unknowns.  To  solve  the  local  saddle  point  problems,  we  use  the  following  Schur 
complement  method: 

1.  Solve  AhUo  =  Whi  in  the  first  step.  For  uniform  meshes.  Ah  is  tridiagonal  for  a 
suitable  ordering  of  the  unknowns,  provided  that  A{x^y)  is  isotropic,  which  is  the 
case  when  A{x^y)  is  diagonal.  Thus,  we  use  direct  solvers  in  this  step. 

2.  Solve  BhAj^^ B^ph  =  Fh  —  BhV-o,  using  a  conjugate  gradient  method.  As  mentioned 
earlier,  even  for  A{x^y)  =  /,  this  matrix  becomes  progressively  ill  conditioned  as 
h  — '  0.  In  the  case  that  A[x,y)  is  piecewise  constant,  the  condition  number  of 
BhAj^^B^  is  at  least  on  the  order  of  the  jump.  We  have  found  that  for  extreme 
values  of  the  jump,  the  conjugate  gradient  method  does  not  always  converge,  even 
if  the  number  of  iterations  is  larger  than  the  size  of  the  system.  We  therefore  use  a 
diagonally  scaled  conjugate  gradient  algorithm,  using  the  diagonal  of  the  matrix  as  a 
preconditioner.  We  remark  that  other  preconditioners  can  also  be  chosen,  such  as  a 
finite  difference  matrix  associated  with  the  discretisation  of  the  elliptic  problem  for 
p.  See  Wheeler  and  Gonzalez  [33]. 

3.  In  the  third  and  final  step,  we  compute  the  solution  of 

AhUi  =  —B^po. 

The  solution  equals  u  =  Uq  -\-  Ui  and  p. 

Since  the  mesh  is  uniform,  the  interpolation  map  from  the  coarse  to  the  fine 
mesh  is  simple.  As  mentioned  in  earlier  sections,  it  is  possible  to  color  the  extended 
subdomains  so  that  the  iteration  in  the  multiplicative  Schwarz  method  is  more  par- 
allelisable.  In  one  of  the  algorithms  we  tested,  we  colored  the  entire  collection  of 
extended  subdomains  using  four  colors.  In  addition  we  have  the  coarse  mesh  prob- 
lem. Figure  2.3  describes  the  four  colors.  We  refer  to  the  resulting  algorithm  as  the 
four  color-five  subspace  multiplicative  Schwarz  algorithm.  Surprisingly,  as  the  results 
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indicate,  this  method  converges  faster  than  a  multiplicative  Schwarz  method  based 
on  the  standard  lexicographic  ordering  of  the  extended  subdomains. 

We  remark  that  in  testing  the  algorithms,  we  found  that  the  inclusion  of  the  coarse 
model  problem  enhanced  the  convergence,  and  provided  a  rate  of  convergence  quite 
independent  of  the  mesh  parameters  n,  n,  and  n,,.  Without  this  coarse  model,  which  is 
quite  an  inexpensive  addition,  the  rate  of  convergence  deteriorated  with  sin  increasing 
number  of  subdomains.  This  is  to  be  expected  for  the  following  reason.  Using  the 
properties  of  the  Green's  function,  we  see  that  the  solution  of  an  elliptic  problem  with 
/(i)  having  support  in  Q,*  C  fi,  can  have  values  significantly  different  from  zero  in  all 
of  Q.  Thus,  in  a  domain  decomposition  method  in  which  there  is  only  a  transfer  of 
information  between  adjacent  subdomains  during  each  iteration,  we  woiild  be  able  to 
devise  a  problem  for  which  at  least  np  iterations  are  required  for  convergence.  Here, 
tIjd  denotes  the  largest  number  of  steps  required  to  transfer  information  between  any 
two  subdomains  in  Q;  see  Widlund  [37].  The  coarse  model  problem  is  an  attempt 
to  include  a  mechanism  to  transfer  information  to  all  the  subdomains,  during  each 
iteration,  i.e.,  a  mechanism  for  a  "global  transportation  of  information". 

For  the  multiplicative  Schwarz  methods,  we  store  the  local  pressures  that  are 
obtained  in  computing  the  projections  onto  the  local  divergence  free  spaces,  for  the 
most  recent  iteration.  However,  we  do  not  use  these  local  pressures  to  compute 
the  residuals.  With  these  stored  local  pressures,  we  are  able  to  compute  the  global 
pressure  ph  once  the  iteration  is  stopped,  in  0{n^)  flops.  Since  the  Schwarz  method 
involves  incrementing  the  current  iterate  by  adding  an  update  obtained  by  solving 
a  subproblem,  the  maginitude  of  the  error  in  the  local  problem  solver  contrctins  the 
accuracy  of  the  Schwarz  method.  In  other  words,  if  the  local  problems  are  accurate 
only  to  within  a  tolerance  of  e,  the  iterates  in  the  Schwarz  method  would,  in  general, 
be  accurate  to  only  the  same  tolerance.  However,  this  need  not  be  the  case  in  the 
acclerated  symmetrised  multiplicative  Schwarz  method,  which  we  did  not  test. 

In  testing  the  algorithms,  we  found  that  for  problems  in  which  the  jump  ratio 
1/J  is  large,  the  relative  error  of  the  computed  pressure  solution  was  larger  thaji  the 
relative  error  of  the  computed  velocity  solution  by  a  factor  on  the  order  of  1/J.  We 
observed  this  phenomenon,  even  when  we  used  direct  methods  to  solve  the  saddle 
point  problems.  We  attribute  this  to  round-off  error. 
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Figure  2.3:  Coloring  a  coarse  mesh  containing  sixteen  subdomains 
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We  remark  that  for  the  additive  Schwarz  method,  we  compute  the  right  hand  side 

9h  =  Y^  PjUh^DF, 
j=Q 

by  adding  the  projections  onto  the  local  divergence  free  spaces,  see  section  2.1.2  and 
section  2.2.2.  Thus,  gh  is  divergence  free,  and  the  iterates  in  the  conjugate  gradient 
algorithm  remain  in  the  divergence  free  space  Ti,  within  the  error  of  the  subproblem 
solvers. 

We  now  present  the  numerical  results  in  table  2.1,  table  2.2  and  table  2.3.  We  list 
the  convergence  factor  for  the  various  mesh  parameters  and  jumps.  The  jumps  listed 
are  the  ratios  1/J;  see  equation  (2.34).  The  convergence  factor  is  the  average  factor  by 
which  the  relative  error  is  reduced,  during  an  iteratioii,  using  the  standard  euclidean 
mesh  norm.  They  are  virtually  the  same  in  the  A  norm  and  in  the  maximum  norm. 
The  exact  solutions  were  randomly  chosen  using  a  random  number  generator  using 
the  uniform  distribution  on  (  —  2,2).  The  number  of  velocity  unknowns  for  n  =  16  is 
480  and  the  number  of  pressure  unknowns  is  256.  For  n  =  24,  there  are  1104  velocity 
unknowns  and  576  pressure  unknowns.  For  n  =  32,  there  are  1984  velocity  unknowns 
and  1024  pressure  unknowns.  For  n  =  40,  there  are  3120  velocity  unknowns  and  1600 
pressure  unknowns.  For  n  =  60,  there  are  7080  velocity  unknowns  and  3600  pressure 
unknowns. 

Conclusions.  The  numerical  results  indicate  that  the  convergence  factor  of  the 
four  color-five  subspace  multiplicative  Schwarz  method  is  about  the  square  of  the 
convergence  factor  of  the  additive  Schwarz  method  with  a  coarse  model.  Thus,  the 
additive  Schwarz  method  with  the  coarse  model  would  require  about  twice  as  many  as 
iterations  as  for  the  four  color-five  subspace  multiplicative  Schwarz  method,  to  reduce 
the  relative  error  by  the  same  factor.  The  results  also  indicate  that  the  lexicographic 
multiplicative  Schwarz  algorithm  with  a  coarse  model  is  slower  than  the  four  color- 
five  subspace  multiplicative  Schwarz  method.  For  all  of  these  methods,  the  resiJts 
indicate  that  the  rate  of  convergence  is  quite  independent  of  the  jump  J,  and  mildly 
dependent  on  the  amount  of  overlap  ratio,  Tionjn.  In  addition,  the  results  indicate 
a  rate  of  convergence  that  is  overall,  quite  independent  of  the  mesh  parameters,  for 
fixed  overlap  ratios. 

Recall  that  our  proofs  on  the  rate  of  convergence,  did  not  make  significant  use  of 
the  coarse  model.  Without  this,  it  is  not  possible  to  prove  that  the  rate  of  convergence 
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Table  2.1:  Results  for  the  additive  Schwarz  method  with  a  coarse  model. 


h-'=n 

Decompositions 

n,  *  n. 

no 

Overlap 
ratio  non,/n 

Jump 

Convergence  factor  ; 

16 

2  *  2 

1 

1/8 

1 

0.3428 

16 

2*  2 

1 

1/8 

10« 

0.4287 

16 

2  *  2 

2 

1/4 

1 

0.2645 

16 

2*  2 

2 

1/4 

10^ 

0.3804 

16 

4*4 

1 

1/4 

1 

0.3434 

16 

4*  4 

1 

1/4 

10« 

0.4057 

16 

4*4 

2 

1/2 

1 

0.3443 

16 

4*  4 

2 

1/2 

10^ 

0.3693 

16 

8  *  8 

1 

1/2 

1 

0.3417 

16 

8*  8 

1 

1/2 

10^ 

0.3494 

24 

4*  4 

1 

1/6 

1 

0.3788 

24 

4*  4 

1 

1/6 

10^ 

0.4106 

24 

4*  4 

2 

1/3 

1 

0.3399 

24 

4*  4 

2 

1/3 

10^ 

0.3405 

24 

8*  8 

1 

1/3 

1 

0.3483 

24 

8*  8 

1 

1/3 

10^ 

0.3589 

32 

4*  4 

1 

1/8 

1 

0.4227 

32 

4*  4 

1 

1/8 

10^ 

0.4436 

32 

4*  4 

2 

1/4 

1 

0.3335 

32 

4*  4 

2 

1/4 

10"^ 

0.3730 

32 

8*  8 

1 

1/4 

1 

0.3528 

32 

8*  8 

1 

1/4 

10^ 

0.3694 

32 

8*  8 

2 

1/2 

1 

0.3492 

32 

8*  8 

2      1/2 

10^ 

0.3520 

40          4*  4 

1 

1/10 

1 

0.4402 

40          4*  4 

1     1/10 

10^ 

0.4843 

40          4*  4 

2 

1/5 

1 

0.3597 

40          4*  4 

2 

1/5 

10^ 

0.3852 

40 

8*  8 

1 

1/5 

1 

0.3721 

40 

8*  8 

1 

1/5 

10^ 

0.3935 

40 

8*  8 

2 

2/5 

1 

0.3325 

40 

8*  8 

2 

2/5 

10^ 

0.3436 

40 

10  *  10 

1 

1/4 

1 

0.3552 

40 

10  *  10 

1 

1/4 

10« 

0.3796 

40 

10  *  10 

2 

1/2 

1 

0.3530 

40 

10  *  10 

2 

1/2 

10^ 

0.3590 

40 

20  *  20 

1 

1/2 

1 

0.3455 
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TaVle  2.2:  Results  for  the  four  color-five  subspace  multiplicative  Schwarz  method. 


\  h~'  =n 

Decompositions 

no 

Overlap 
ratio  nonjn 

Jump 

Convergence  factor 

16 

2  *  2 

1 

1/8 

1 

0.1652 

16 

2*  2 

1 

1/8 

10^ 

0.1627 

16 

2*2 

2 

1/4 

1 

0.0548 

16 

2*  2 

2 

1/4 

10^ 

0.0538 

16 

4  *  4 

1 

1/4 

1 

0.0596 

16 

4*  4 

1 

1/4 

10^ 

0.0612 

16 

4  *  4 

2 

1/2 

1 

0.0276 

16 

4*  4 

2 

1/2 

10^ 

0.0335 

16 

8  *  8 

1 

1/2 

1 

0.0451 

16 

8*  8 

1 

1/2 

10« 

0.0543 

24 

4*  4 

1 

1/6 

1 

0.0939 

24 

4*  4 

1 

1/6 

10« 

0.0987 

24 

4*  4 

2 

1/3 

1 

0.0500 

24 

4^  4 

2 

1/3 

W 

0.0514 

24 

8*  8 

1 

1/3 

1 

0.0530 

24 

8*  8 

1 

1/3 

10^ 

0.0535 

32 

4*  4 

1 

1/8 

1 

0.1584 

32 

4*  4 

1 

1/8 

10^ 

0.1629 

32 

4*  4 

2 

1/4 

1 

0.0576 

32 

4*  4 

2 

1/4 

10^ 

0.0576 

32 

8*  8 

1 

1/4 

1 

0.0622 

32 

8*  8 

1 

1/4 

10^ 

0.0621 

32 

8*  8 

2 

1/2 

1 

0.0360 

32 

8*  8 

2 

1/2 

10« 

0.0402 

40 

4*  4 

1 

1/10 

1 

0.2535 

40 

4*  4 

1 

1/10 

10« 

0.2552 

40 

4*  4 

2 

1/5 

1 

0.0723 

40 

4*  4 

2 

1/5 

10^ 

0.0688 

40 

8*  8 

1 

1/5 

1 

0.0748 

40 

8*  8 

1 

1/5 

10^ 

0.0748 

40 

8*  8 

2 

2/5 

1 

0.0406 

40 

8*  8 

2 

2/5 

10^ 

0.0391 

40 

10  *  10 

1 

1/4 

1 

0.0644 

40 

10  *  10 

1 

1/4 

10« 

0.0641 

40 

10  *  10 

2 

1/2 

1 

0.0360 

40 

10  *  10 

2 

1/2 

10^ 

0.0392 

40 

20  *  20 

1 

1/2 

1 

0.0748 

60 

10  *  10 

1 

1/6 

10<^ 

0.0884 

60 

10  *  10 

2 

1/3 

1 

0.0683 
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Table  2.3:  Results  for  the  lexicographic  multiplicative  Schwarz  method. 


h-'   =n 

Decompositions 

Tlo 

Overlap 
ratio  UoTiJu 

Jump 

Convergence  factor 

16 

2*2 

1 

1/8 

1 

0.1613 

16 

2*  2 

1 

1/8 

10^ 

0.1466 

16 

2*2 

2 

1/4 

1 

0.0558 

16 

2*  2 

2 

1/4 

10^ 

0.0486 

16 

4  *  4 

1 

1/4 

1 

0.0836 

16 

4*  4 

1 

1/4 

10^ 

0.0888 

16 

4*4 

2 

1/2 

1 

0.0545 

16 

4*  4 

2 

1/2 

w 

0.0631 

16 

8  *  8 

1 

1/2 

1 

0.1143 

16 

8*  8 

1 

1/2 

10« 

0.1253 

24 

4*  4 

1 

1/6 

1 

0.0926 

24 

4*  4 

1 

1/6 

10^ 

0.1035 

24 

4*  4 

2 

1/3 

1 

0.0635 

24 

4*  4 

2 

1/3 

10« 

0.0772 

24 

8*  8 

1 

1/3 

1 

0.1048 

24 

8*  8 

1 

1/3     10« 

0.1236 

32 

4*  4 

1 

1/8 

1 

0.1550 

32 

4*  4 

1 

1/8 

10^ 

0.1727 

32 

4*  4 

2 

1/4 

1 

0.0773 

32 

4*  4 

2 

1/4 

10^ 

0.0911 

32 

8*  8 

1 

1/4 

1 

0.1027 

32 

8*  8 

1 

1/4 

10^ 

0.1170 

32 

8*  8 

2 

1/2 

1 

0.1132 

32 

8*  8 

2 

1/2 

10^ 

0.1245 

40 

4*  4 

1 

1/10 

1 

0.2499 

40 

4*  4 

1 

1/10 

10*^ 

0.2644 

40 

4*  4 

2 

1/5 

1 

0.0735 

40 

4*  4 

2 

1/5 

10^ 

0.0985 

40 

8*  8 

1 

1/5 

1 

0.0944 

40 

8*  8 

1 

1/5 

10^ 

0.1013 

40 

8*  8 

2 

2/5 

1 

0.0987 

40 

8*  8 

2 

2/5 

10^ 

0.1137 

40 

10  *  10 

1 

1/4 

1 

0.1011 

40 

10  *  10 

1 

1/4 

w 

0.1120 

40 

10  *  10 

2 

1/2 

1     ' 

0.1159 

40 

10  *  10 

2 

1/2 

w 

0.1299 
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is  independent  of  the  coarse  mesh  parameter  H  =  l/n, ,  since  without  a  mechanism  for 
a  "global  transportation  of  information"  during  each  iteration,  the  rate  of  convergence 
of  the  algorithms  will  depend  on  H.  We  have,  so  far,  been  unable  to  prove  that  the 
rate  of  convergence  is  independent  of  H.  We  mention  that,  for  standard  finite  element 
discretisations  of  elliptic  problems,  Dryja  and  Widlund  [13]  have  proved  that  the  rate 
of  convergence  of  the  additive  Schwarz  method  is  independent  of  both  H  and  h,  in 
both  two  and  three  dimensions. 
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Chapter  3 

Iterative  refinement  methods  for 
Raviart-Thomas  elements. 


In  this  Chapter,  we  discuss  iterative  methods  to  solve  the  discrete  linear  systems 
arising  from  the  discretisation  of  mixed  formulations  of  elliptic  problems  using  the 
Raviart-Thomas  spaces  on  a  repeatedly  refined  mesh.  These  iterative  methods  are 
referred  to  as  the  FAC  and  AFAC  algorithms,  in  a  similar  context  for  the  discreti- 
sation of  standard  elliptic  problems.  First,  we  describe  the  iteratively  refined  mesh, 
and  the  Raviart-Thomas  spaces.  Following  that  we  describe  the  FAC  and  AFAC  al- 
gorithms and  present  some  theoretical  bounds  for  their  rates  of  convergence.  We  also 
include  some  quantitative  bounds  for  the  rate  of  convergence  of  the  FAC  algorithm, 
in  the  special  case  of  discretisation  using  the  lowest  order  Raviart-Thomas  spaces  on 
triangles. 

3.1      Definition  of  the  refined  meshes  and  spaces. 

Let  ^0  =  Q  C  R^  he  a.  polygonal  domain  triangulated  by  a  uniformly  shape  regular 
mesh  of  width  ho-,  denoted  by  r^°{Q.o).  Let  Qi  C  fio  be  a  polygonal  subregion  of  the 
triangulation  where  we  wish  to  refine  the  mesh.  If  there  is  neccesity  for  still  further 
refinement,  say,  another  /  —  1  levels  of  refinement,  we  introduce  a  sequence  of  nested 
polygoncd  regions  where  we  intend  to  refine  the  mesh,  repeatedly: 

Before  we  define  the  repeatedly  refined  mesh,  we  introduce  a  standard  multilevel  grid 
defined  on  Qq.  We  consider  a  sequence  of  refined  uniformly  shape  regular  meshes,  of 
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mesh  widths, 

hi  <  /i;_i  <  ■  ■  ■  <  hi  <  ho, 

with  the  corresponding  partitions 

t'^I^o)  C  T'"(fio)  C  •  •  •  C  r'"-'(no)  C  T'"(f]o).  (3.1) 

We  assume  that  the  vertices  of  the  elements  in  the  triangulations  (3.1)  never  lie  on 
the  interior  of  an  edge  of  another  element.  By  a  refinement  we  mean  the  following: 
Each  element  if  of  a  coarser  mesh  is  divide  into  two  or  more  smaller  elements,  such 
that  the  vertices  of  the  smaller  elements  do  not  lie  on  the  interior  of  an  edge  of  a 
neighbouring  element  in  the  refined  mesh.  In  a  standard  case  of  refinement,  as  in  the 
case  of  many  multigrid  methods,  each  element  in  the  coarser  mesh  is  divided  into  four 
subelements  by  connecting  the  midpoints  of  the  coarse  mesh  element.  In  this  case, 
the  mesh  width  of  the  refined  mesh  is  reduced  by  a  factor  of  2. 

For  i  —  1,...,/,  the  uniformly  shape  regular  mesh  t^'{Q,q),  of  width  hi  <  hi-i, 
is  assumed  to  be  constructed  by  refinement  of  the  uniformly  shape  regular  mesh 
t''"-'(J7o)-  We  assume  that  fi,,  for  i  =  1,...,/,  can  be  expressed  as  a  union  of  ele- 
ments in  t'''-'(Ho)-  To  obtain  the  iterattvely  refined  mesh,  we  consider  restrictions  of 
each  mesh  t^'{Q.o)  to  the  subregion  Q,  to  obtain  r'''(n,). 

DEFINITION.     We  define  the  repeatedly  refined  mesh  t^° '''(f^)  as  consisting  of  those 

elements  which  are  in  any  of  the  various  refined  meshes  T'''(r2o)  restricted  to  fi^  — fii+i, 
where  we  define  fi;+i  =  0.  i.e., 

> '"(fi)  =  uLoT'"'(a-o,+i). 

REMARK.  In  the  case  of  the  repeatedly  refined  mesh,  the  vertices  of  elements  in  the 
refined  regions  such  as  along  dQ,i,  may  lie  on  the  interior  of  edges  of  coarser  elements 
in  neighbouring  regions.  See  Figure  3.1.  Also,  the  regions,  Qi  need  not  be  connected, 
and  this  can  lead  to  increased  parallelism  in  the  numerical  algorithms  to  be  studied. 

DEFINITION.     We  define  the  refined  mesh  Raviart-Thomas  spaces: 

Xk,....m{^)  =  ^;,„(f^o)  +  X;,,  (no  +  •  •  ■  +  j^h,_,  (a-:)  +  J^^,(a), 
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Figure  3.1:  A  repeatedly  refined  rectangular  mesh  with  2  levels  of  refinement. 
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f)i 

Qo 

Oj 

fil 

fio 

where  Xhi{^i)  are  the  Raviart-Thomas  velocity  spaces  with  zero  normal  trace  on  5f2,-, 
as  used  in  the  discretisation  of  the  Neumann  problem  on  Q,.  Of  course, 


Xho h,{^)  ^  Ho,an{div,Q). 


Similarly,  we  define 


where  y|,,(fii)  are  the  Raviart-Thomas  pressure  spaces  with  zero  mean  value  on  Cli  as 

used  in  the  discretisation  of  the  Neumann  problem  on  Q,.  Thus  Yh,^ h,{^)  is  also  a 

subset  o{  L^{n)/R. 

We  use  the  refined  mesh  Raviart-Thomas  spaces  {Xho,..,hii^)jYho ;i,(f^)))  de- 
fined on  the  repeatedly  refined  mesh,  to  obtain  our  discretisation  of  saddle  point 
problem  (1.7).  Thus,  our  discrete  problem  is: 


Find  Uh  e  A'^^,...,h,(fi),  ph  G  Yf^_,„^h,i^)  such  that 
a{uK,v)  +  b{v,Ph)    =    0,  W  e  Xh,,...,hX^) 

.    t(^^,g)  =   Snf{^M^)dx,  VqeYho.....h,{^), 


(3.2) 


where  a(., .)  and  6(., .)  are  defined  by  equation  (1.12).  For  well  posedness,  we  need  to 
check  the  following  two  conditions. 
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1.  First,  we  need  to  check  that  a{., .)  is  positive  definite  on 

Xho />,;o(n)  =  {UH  e  Xh, h,{^)  :  b{uH,q)  =  0,     Vg  e  n,....,;.,(f2)}, 

equipped  with  the  H{div,Q)  norm,  with  a  constant  independent  of  h.  We  will  show 
that 

if  Uh  e  Xho hr.oi^)  tten  V  •  ■u;,  =  0, 

i.e.,  discrete  divergence  free  implies  divergence  free  in  the  i^(f2)  sense.   Then,  a{.,.) 
will  be  positive  definite  on  Xha,...,hi;o{^)  with  constant  a  since 

Xho Hrfii^)CHo,anid^v°,^)^ 

and  since 

aiv,v)  >  a||t7||^^,,^^^^,       WeH{div',n), 

see  Chapter  1,  equation  (1.13). 

So  consider,  Uh  G  A';,;,  ..,/,, (fi).  Then,  by  definition, 

Hh  =  Uo  + ■  ■  •  +  ui,       where    u^  G  XH,.(fii),Vi, 

and 

V  •  •iT;,  =  V  ■  Wo -f hV-u;,       where    V  •  Ui  G  1^,(^0, Vi. 

By  definition  of  Yh^ hi{^)  ^^  the  sum  of  the  Raviart-Thomas  pressure  spaces  defined 

on  the  refined  meshes,  it  follows  that  V  •  u/,  G  Yh^^,„^h,{^)-  This  implies  that 

i{b{uh,q)  =  0,    VqeYh, ;.,(n)then    f  \V  ■  Unl^dx  =  0. 

Jn 

Thus  discrete  divergence  free  implies  divergence  free  in  the  L^  sense. 

2.  Next,  we  need  to  check  that  b{.,.)  satisfies  a  uniform  inf-sup  condition  on  the  re- 
fined mesh  Raviart-Thomas  spaces  {Xho_,„^h,{^),Yhf, h,(fi))-  Given  qn  G  lh^,...,/,,(n), 

we  will  show  that  there  exists  a  u /,  G  Xho....,h,{^)j  satisfying  V  •  Uh  =  qh,  and 

ll^''llH(<f.-v,n)  <  c(no,...,f^i)lk'»IU'(n)'  (3-3) 

for  some  positive  constant  c{Clo,  ■  ■  ■  ,^i)  independent  oiho,. . .  ,hi.  From  this  it  follows 

that 

b{wh,qh)  .  b{uH,qh) 


sup  n > 


^^ex,, ^,(n)-{o}  lPhllH(div,n)ll9''|U=(n)        lPh||^(di„,n)ll9MlL^(n) 
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_       /n(V  •  uh)qhdx       _  llg^llL^(n)  y  1 

Il'"'>ll^(d.-v,n)ll9''ll^'(n)        \\^h\\H{d^v,u)\MLHO)  ~  c(51o,  •  •  ■  ,f^;)' 
which  is  the  uniform  inf-sup  condition  with  constant  l/c{Qo,. . .  ,Qi). 

Now  we  prove  condition  (3.3).    Let  P,i7  denote  the  L^  projection  onto  Yh,{Cli), 
for  1  =  0 ,/.  We  will  find  an  L^{^)  orthogonal  decomposition  of 

Qh  =  qo  +  •  •  •  +  qi, 
with  Qi  G  yji, (Q,).  For  i  =  0,  we  let 

go  =  Po.L^qh, 

and  for  z  =  1, ...,/,  we  let 

Note  that,  for  1  >  z, 

qn-Y.  9.  =  (^  -  ^0,1^  -  E  P:.l^PUl^  ■  ■  ■  Po,L^)qk 

3=0  J=l 

=  -Pi-i,L'  ■  ■  ■  Po'.L^qh, 

and  thus, 

i-l 

5.  =  P^.L^iqh  -  2Zgj)- 

j  =  0 

Since  both  qh  and  go  have  mean  value  zero,  it  follows  that  go  is  the  L^  projection  onto 
yho{^o)®'Po{^o)-  By  discontinuity  of  the  pressure  spaces  across  element  boundaries, 
and  the  fact  that  Yh^{Qo  —  Qi)  =  y/,„,...^/i,(no  —  f^i),  it  follows  that 

qh  -  qo  =  0,        on  fio  -  f^i- 

Since,  the  characteristic  function  of  fii,  denoted  by  Xfii  >  is  in  y/ij,(no)  @  PoC^o))  it 
follows  that 

Qh  —  qo  has  mean  value  zero  on  Vli. 

Similarly,  for  i  =  2, . . . ,  /,  we  obtain  by  induction,  that 

t-i 
<?/.  -  X^  g^  =  0       on  Qo  -  Q,, 

j=0 
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and 

t-i 


Thus  Qi,  which  is  the  L^  projection  of  g/,  —  Yl)=oQj  o°^°  ^.(^>'))  ^^  satisfy 

i 

Qh  -Yll^  -  ^       on  Jlo  -  f^i+i, 

j  =  0 

by  using  the  discontinuity  of  the  discrete  pressures  across  element  boundaries,  and 
that 

n.(a  -  a+i)  =  n„ hM  -  a+i). 

Note  that, 

/-I  /-I 

j=o  j=o 

Thus, 

g/.  =  go  + \-  qi- 

We  show  that  {g,}  are  mutually  orthogonal  as  follows.    First,  for  i  =  1,  we  will 
show  that  (Jo  ±  gj.  By  construction, 

qH-qoeY,,ino)^,  (3.4) 

and  thus  qh  —  qo  has  support   in  fli.  Since  gi  is  the  L^  projection  of  q^  —  go  onto 
F^^(Qi),  it  follows  that 

I  {qh  -  qo  -  q  -  l)qodx  =  [   {qh  -  qo  -  q  -  l)qodx  =  0,  (3.5) 

because  golfii  ^  i'/ii(^i)-  By  using  equation  (3.4)  and  equation  (3.5),  we  obtain  that 

/     qoqydx  =  0. 
Jo, 

Next,  we  show  that  if  1  <  z,  then  g,+i  is  orthogonal  to  qo,...,qi.  By  definition, 

g,+i  is  the  L^  projection  of  qn  -  E}=o  9j'  which  has  support  in  fi.+i,  onto  y;^.^l(^i+l). 

Thus, 

i+i  .  i+i 

/    (g^^  -  I]  gj)gfc^2:  =    /       {qh-Yl  9:)qkdx  =  0, 
for  k  <  i,  since  g^ln.+i  e  yh^(nj),  for  j  >  k.  We  decompose 

/      (?/>  -  H  gjJg'c^^  =  /      {qh-Yl  (iMkdx  +  /      (-g.+i)gfcc^i. 
•/n.+i  j=o  •'"•+■  i=o  "'"•+■ 
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The  first  part  is  zero  by  the  definition  of  g,,  since  qk\o,^i    £  ^'^/i, (^t)-  Therefore,  the 
second  part  satisfies 

0  =  /       {-qr+i)qkdx, 
which  shows  that  q^  -L  q,+i  for  k  <  i.  Thus, 

qh  -  qo  + \-  qi,        where  g,-  G  Yh,{Cli), 

is  an  orthogonal  decomposition  of  qn- 

Since  the  triangulation  T^'{Cli)  is  uniformly  shape  regular  on  fi,,  we  know  that  the 
uniform  inf-sup  condition  holds  on  the  local  Raviart-Thomas  spaces  for  the  Neumann 
problem,  (A'/,, (Q,),  y/,^(fi,)).  Thus,  there  exists  a  positive  constant,  Ci{Q,)  such  that 
given  q,  G  }/,,(n,),  there  exists  r,  G  A'"/,,(n,)  such  that  V  •  v^  —  q,  and 


By  defining 
we  obtain, 
and 


V  -4  =  9o  + r  qi  =  qh, 


:';=0  i=0 


<  max(co,...,c;)\//  -f  l\\qh\\mQ)- 

Note  that,  the  constants  Co,...,c/  are  independent  of  ho,....,hi.  Therefore,  the  re- 
peatedly refined  Raviart-Thomas  spaces  satisfy  an  inf-sup  condition,  with  a  constant 
independent  of  /iq,  • . . ,  /z;,  and  possibly  depending  on  /,  the  number  of  levels  of  refine- 
ment, and  the  geometry  of  the  refined  regions. 

Since  the  repeatedly  refined,  composite,  spaces  lead  to  stable  discretisations,  the 
error  in  discretisation  is  given  by 

where  {u,p)  is  the  solution  of  the  continuous  problem  and  (u;,,p/,)  is  the  solution  of 
the  discrete  problem.  We  assume  that  our  choice  of  refinement  leads  to  a  small  error. 
Throughout  the  rest  of  this  Chapter,  we  consider  iterative  methods  to  solve  saddle 
point  linear  system  (3.2). 
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3.2      FAC  and  AFAC  methods  for  Raviart- Thomas 
elements. 

The  fast  adaptive  composite  grid  methods  (FAC),  and  the  asynchronous  fast  adap- 
tive composite  grid  methods  (AFAC)  were  originally  introduced  as  iterative  methods 
to  solve  positive  definite  linear  systems  arising  from  the  standard  discretisations  of 
elliptic  problems  using  repeatedly  refined  meshes;  cf.  McCormick  and  Thomas  [26], 
Mandel  and  McCormick  [24],  Widlund  [36]  and  Dryja  and  Widlund  [14].  The  two 
level  case  is  discussed  in  Mandel  and  McCormick  [24]  and  also  Bramble,  Ewing,  Pas- 
ciak  and  Schatz  [6].  We  refer  to  Widlund  [36]  and  Dryja  and  Widlund  [14]  for  the 
analysis  of  the  many  level  case.  The  FAC  and  AFAC  algorithms  are  actually  the 
multiplictive  and  additive  versions  of  the  Schwarz  method,  respectively,  applied  to 
special  subspaces. 

As  discussed  in  chapter  2,  our  application  to  the  saddle  point  problem  is  accom- 
plished in  three  steps.  The  second  step  is  where  the  Schwarz  iterative  methods  are 
applied,  and  this  step  is,  computationally,  the  most  expensive  step.  Steps  1  and  3,  are 
virtually  the  same  for  the  FAC  and  AFAC  methods,  although  less  work  is  involved  in 
the  third  step  of  the  FAC  method.  The  steps  are: 

1.  Find  Uh,j  €  Xho,...,h,{^)  such  that 

BhUh.i  -  Fh. 

2.  Find  Uh,DF  ^  ^h(,....,hi  such  that 

a{uh,DF,v)  =  <i{-Uh.i,v),       We  Tiho hn 

where  Ho i  =  Xh, h,{^)  H  H{div°,n). 

3.  Determine  the  pressure  ph- 

3.2.1      Step  1:   Reduction  to  a  divergence  free  problem. 

To  find  a  discrete  velocity  Uhj  with  the  same  divergence  as  the  discrete  velocity 
solution  Uh,  we  sequentially  solve  problems  on  each  of  the  nested  regions  determining 
a  normal  trace  compatible  with   f{x)  on  the  next  subregion.     First,  we  solve  the 
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following  problem  on  the  coarsest  mesh  t^°{CIo): 

Find  Uho,i  e  Xhoi^o),  Pho.i  ^  ^^0(^0)  such  that 
<      a{uHo,i,v)  +  b{v,pf,^j)     =     0  WeXhoi^o)  (3.6) 

,     b{u,„i,q)  =     Jafix)qdx       Vg  G  Fh„(no) 

Substituting  q  G  y/n,(r2o)  with  support  in  f^o  —  Jli,  we  see  that  the  solution  of  equa- 
tion (3.6),  Uh^j  satisfies 

"^  •  Uho.i  -  "^  •  Uh,        on  ^0-^1, 

where  Uh  Js  the  discrete  velocity  solution,  problem  (3.2)  (repeatedly  refined  saddle 
point  problem).  Note  that,  on  fij,  the  divergence  of  Uh  —  'W/io,/  will  not  be  zero  in 
general.  However,  Uho,l  will  have  normal  trace  on  dQi  compatible  with  the  divergence 
constraint,  i.e., 

/     ^■Uk^,jdx=   /     f{x)dx=   /      ')„UhaidS:,. 

Thus,  f{x)  —  V  •  Uho.l  has  mean  value  zero  on  Qj,  and  so  zero  normal  trace  on  dQ.i 
is  compatible  with  it.  Thus,  we  could  pose  a  problem  on  Qi  using  the  new  residual, 
and  repeat  the  procedure,  sequentially,  on  ni,...,Q;.  We  define 

I 

j  =  0 

where  each  Uh,j  is  obtained  by  solving  the  following  local  problem,  for  i  =  !,...,/: 

Find  Uh,.i  e  XhS^x),  Ph,,i  G  i^.l^^,)  such  that 

a{uH,j,v)^b{v,Ph,j)     =     0  \/veXH,{Q,)       (3.7) 

b{uH.j,q)  =     Ja[f{x)-ZTJo'^-UH,.i]qdx      Vg  G  ^/..(a) 

Note  that,  by  construction,  for  i  =  0, ...,/  —  1, 

/    Xn.+,(V--u^.,/)<fx  =  /    Xu,^Af{x)-Y.y  ■Uh,,i]dx 

=     /  Inilh.jdS:,. 

Thus,  f{x)  -  E}=o  ^  ■  ■"/>;,/  h^s  mean  value  zero  on  fii+i,  which  makes  problem  (3.7) 
for  I  +  1  well  posed.  Note,  that,  for  i  =  0, . . . ,  /, 

j'-i 
V  -Uh,,!  =  V-Uh-^y-Uh,,i,        on  Cli  -  Q.+i, 

j  =  0 
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where,  Cli+i  =  0.  Thus,  we  obtain 

i.e.,  BhUhj  =  Fh.  The  equation  satisfied  by  {uh  —  ■Uh,i,Ph)  is: 

f  a{uh,DF,v)  +  b{v,Pk)     =    0-a{uh,i,v)     W  e  Xh^,...^h,{^) 
I    b{uH,DF,q)  =     0,  Vgen, H,{^), 


(3.8) 


where  u/,  ci?  =  Uh  —  Uh,i-  Thus,  Uh,DF  is  divergence  free  and  is  the  solution  of  the 
following  equivalent  problem: 

J  Find  Uh.DF  €  'Hho,...,h,  such  that  ,„  .. 

\   a{'Uh.DF,v)  =  -a{uh,i,v)        \fv  e 'Hko....,h,, 

where  ?^,, ,,  =  A',„ ,,(Q)  n  H{div°,n). 

3.2.2      Step   2:     Solution   of  the   divergence   free   symmetric 
positive  definite  problem. 

We  solve  problem  (3.9)  by  the  use  of  the  multiplicative  or  additive  Schwarz  methods, 
or  as  referred  to  in  this  context,  the  FAC  and  AFAC  algorithms,  respectively.  Once 
'^h,DF  is  obtained,  the  discrete  velocity  solution  Uh  is  determined  by  letting 

First,  we  introduce  some  notation.  For  z  =  0, .  .  .  ,  /,  let 

m,  =A\(Q.)n^o,an.(<iii'°,a). 

These  are  divergence  free  subspaces  of  the  Ravi  art -Thomas  velocity  spaces  defined 
on  the  refined  meshes,  T^°{no),  ■  ■  ■  iT^'i^i)-  The  functions  in  these  subspaces,  when 
extended  by  zero  outside  Q,  to  the  whole  domain  Q,  lie  in  the  global  space  Ti/i^, ...,;,,. 
These  are  the  subspaces  used  in  the  formulation  of  the  FAC  and  AFAC  methods. 

Basis  for  the  divergence  free  subspaces  'HhQ,...,hii'Hhoi  •  ■  •  t'^hi  aje  not  available. 
To  implement  the  Schwarz  methods,  we  need  to  compute  the  projections  onto  the 
divergence  free  subspaces.  We  compute  these  by  solving  saddle  point  subproblems. 
For  i  —  0,...,/,  let  P,  denote  the  orthogonal  projection  in  the  a(.,.)  inner  product 
onto  Tih,-   To  obtain  each  projection,  we  solve  a  saddle  point  subproblem  using  the 
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basis  for  the  appropriate  local  refined  mesh  Raviart-Thomas  spaces.  For  instance,  to 
find  PiVJ,  for  2  =  0, ...  ,  /,  we  solve: 

Find  Wh,  G  Xh,{^')j  ^h,  G  y^, (fJ,)  such  that 

a{wh,,v)  +  b{v,Sh,)     =     a{w,v)    ViT  G  A'/..(Q.)  (3.10) 

b{wh„q)  =     0  Vgen,(a) 

Thus  PiW  =  Wh,  G  Tih,  ■ 

We  now  state  and  prove  a  result  concerning  the  existence  of  a  partition  for  the 
divergence  free  Raviart-Thomas  subspaces.  This  result  will  be  used  to  prove  conver- 
gence results  for  the  FAC  and  AFAC  methods  applied  to  elliptic  problems  discretised 
by  the  Raviart-Thomas  mixed  finite  element  method. 

Lemma  22  Let  Tih^,  ■  •  ■  ■,'Hh,  be  divergence  free  subspaces  'of'Hho,...,h,i  defined  on  the 
refined  meshes,  as  described  before.  Then,  there  exists  a  positive  constant  fx,  inde- 
pendent of  ho, . . .  ,hi  such  that  for  every  iJ^  G  'HhQ,...,hi,  there  exists  a  partition  with 
^1  G  Tih, ,  for  1  =  0,...,/,  satisfying 


ind 


Vh  =  Vo  +  Vi  ~  ■  ■  ■  +  Vi, 


I 
^a{^^.,^^I)  <  ^^a[vh,Vh). 

1  =  0 


PROOF   OF   LEMMA.     Given  Uh  G  Tih, h,,  we  define  iT,  for  i'  =  0, ...,/-  1  by: 

-  ^  i  ^h  ~  E}=o  ^j      on  Q,  -  fi,+i 
""'       \  E^-g,  on  Q,  +  i 

where  g^  =  fn{uh  —  ^}=o^j)  on  dQ,+i,  and  E^'  denotes  the  Raviart-Thomas  diver- 
gence free  velocity  extension  onto  Vh,{^j-fi)  r\  H {div° ,  ^,+1),  as  given  by  the  extension 
theorem  in  chapter  1.  We  define  u/  as: 

;-i 
ui  =  Uh  —  2_j  '^ji        on  r2(. 

J  =  0 

Note  that,  by  the  normal  trace  lemma 

t-i 

1=0 
and  by  the  extension  theorem 

ll^'''5.ll//(<i,.,n,+,)  <  c,-||5,||;^-,/2(sn.+,)- 
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Combining  these,  we  obtain, 

»-i  «-i 

j=Q  j=0 

for  some  positive  constant  c,-  independent  of  the  mesh  parameter  hi. 
For  i  =  0,  we  obtain 

a{uo,Uo)  <  cla{uh,Uh), 

for  some  positive  constant  cq  independent  of  ho-  By  induction,  we  obtain,  for  i  = 
1,...,/: 

a{u„u,)  <  (1  +  cl){l  +  2cl)  •••(!  +  icf_i)(i  +  iy,an,{uh,Uh). 

Summing  over  all  z,  we  obtain  an  estimate  for  /x': 

^'  <  [{1  ^  (I  ^  l)cf  }{1  4-  Icl,}  .•.{!  +  2c^}{l  +  4}  -  1].      D 

REMARK .  The  preceeding  result  implies  that  for  a  fixed  number  of  levels  /,  the  con- 
stant fi  is  independent  of  the  mesh  parameters  ho, . . .  ,hi.  Since  we  used  the  extension 
theorem  to  obtain  a  bound  for  ji,  it  could  possibly  depend  on  the  subdomain  geom- 
etry. However,  we  will  present  a  result  later  on  that  gives  a  bound  for  the  rate  of 
convergence  of  the  FAC  method  which  is  independent  of  the  geometry  of  the  subdo- 
mains  f^,. 
We  consider  applications  of  this  lemma. 

Theorem  5    The  error  propagation  map  of  the  FAC  algorithm,  i.e.,  the  multiplicative 

Schwarz  method  using  the  subspaces  Tiho,'Hh^, . . .  ^Tih,  of  Hh^ hi,  has  a  convergence 

factor  p  which  is  independent  of  ho, . . .  ,hi. 

PROOF.  We  use  Lemma  19  and  Lemma  22  to  obtain  the  result.  Later  on  we  present 
another  proof  of  this  result.  □ 

The  following  result  concerns  a  variant  of  the  AFAC  algorithm,  cf.  Dryja  and  Widlund 
[14]  and  Widlund  [36].  The  standard  version  of  the  AFAC  algorithm  involves  the  use 
of  projections  P-^-^,  which  are  projections  onto  the  Wj+i  PI  Td.  However,  we  consider 
only  a  variant  of  that  algorithm,  which  involves  Pi,  the  projections  onto  Tii. 

Theorem  6  Let  P,  denote  the  a(.,.)  orthogonal  projection  onto  Tih-.   Then, 

P  =   Po  +  Py  +  ---  +  Pl, 
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is  symmetric  and  positive  definite  in  the  a(.,.)  inner  product.  Furthermore,  there 
exists  positive  constants  61,62,  independent  of  ho, .  ■  ■  ihi,  such  that 

81  a{uh,Uh)  <  a{Puh,Uh)  <  ^2  a{i!h,Uh)- 

PROOF  OF  UPPER  BOUND.  P  is  symmetric  since  it  is  a  sum  of  orthogonal  projections, 
each  of  which  is  symmetric,  in  the  a{.,.)  inner  product.  An  upper  bound  for  ^2  is  /-f  1, 
since  the  norm  of  each  projection  is  bounded  by  one.  We  can  show  that  this  upper 
bound  cannot  be  improved  as  follows.  Let  Uh^  G  Tin,  H  Tiha,  which  will  be  nonempty 
as  ho  is  made  smaller.  Then,  PiUho  ~  ^Ki  for  each  i,  and  so  the  upper  bound  is  /+  1. 
D 

PROOF  OF  LOWER  BOUND.  To  obtain  the  lower  bound,  we  use  Lemma  15,  (Lions' 
lemma),  and  Lemma  22  to  obtain  a  constant  n  independent  of  ho,  •  •  •  ,hi  such  that: 

—  a{v.v)  <  a{Pv,v).        □ 

Next,  we  prove  a  lemma  estimating  the  condition  number  of  the  symmetrised F AC 
method.  Following  that  we  give  an  application  using  the  strengthened  Cauchy  inequal- 
ity to  give  a  bound  for  the  condition  number  which  is  independent  of  the  geometry  of 
the  subregions  fi,,  and  the  mesh  parameters  ho,..., hi.  See  Dryja  and  Widlund  [36], 
and  Mandel  and  McCormick  [24]. 

DEFINITION.     By  the  symmetnsedFAC  algorithm,  we  mean  the  multiplicative  Schwarz 
method  applied  to  solve  equation  (3.9),  using  defect  correction  on  the  subspaces 

in  that  order. 

In  the  following,  we  describe  the  steps  in  more  detail  and  introduce  some  no- 
tation. Let  a{u,v)  denote  the  bilinear  form  associated  with  the  symmetrised  TAC 
preconditioner,  and  let  /(.)  denote  a  linear  functional.  Then,  the  solution  u  of 

a{il,v)  =  f{v),       \/venho....,hn 

is  given  by 

u  =  ui-i-  iii^i  -^  ■  ■  ■  +  Ui  +  uo  +  ul  +  •  ■  ■  +  ii7_i  +  u^,  (3.11) 

where 
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•  ui  £  Tihi  is  the  solution  of 

a{ui,v)  =  f{v),      W  G  Tihi; 

•  ui-i  G  'Hhi_i  is  ^^^  solution  of 

a{ui  +  ut-uv)  =  f{v),      VvG7^h,_,; 


•  Til  G  "W/,!  is  the  solution  of 

a{ui  -t  ui-i  + h  ui ,  z7)  =  /({T),       W  G  ?^hi ; 

•  uo  G  TYho  is  the  solution  of 

a{ui  -f  J;_i  + ^  Ui  +  Uo,v)  =  f{v),       W  G  T^^o; 

•  tx^  G  "H^i  is  the  solution  of 

a{ui  -  -u/-:  H 'r  iii  +  Uq  +  ul.v)  =  f{v),       Vv  G  Ti/n; 


•i7;*_j  G  "W/i,.]  is  the  solution  of 

a{ui  -r  tli_i  -r rUj  -rtro  +  -u^  +  ---  +  ^-ijt^)  =  f{v),       Vi7  G  "W/,,,! ; 


anc 


•  u*  G  T^h,  is  the  solution  of 

a{ui  +  tr;_i  + \-  ui  +  uo  +  ul  + h  iL^_y  +u^,v)  =  f{v),      VvG  'W^,. 

The  following  Lemma  describes  another  representation  for  the  bilinear  form  a{u,v) 
associated  with  the  symmetrised  FAC  algorithm.  See  also  Dryja  and  Widlund  [36] 
and  Mandel  and  McCormick  [24]. 
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Lemma  23    The  bilinear  form  d{u,  v)  has  the  following  representation: 

I 

a(^r,i;)  =  X^M.(u,i;),  (3.12) 

t  =  0 

where 

Mi{u,v)  =  a{Piii,Piv) 

M,{u,v)  =  a(P,//^,^^u,P,if^.^,{r)      /orz  =  0,...,/-  1, 

Pi  denotes  the  a(.,.)  orthogonal  projection  onto  7ih,,  o.nd 

Hq      u  is  defined  for  u  G  'Hho,..,h,  (^nd  0  <  z  <  /  —  1  hy: 

H'^^^^u  e  Tih,  +  ■  ■  ■ -^ n^,, 

H'^^^^u{x)  =  u{x),  /or  X  e  fio  -  n,+i, 

I.e.,  Hq      u  is  the  same  as  u   everywhere  except  on  ^,+1    where  it  is  the  harmonic 
extension  of  u,  in  the  sense  of  Tih^  n?i/,,^j,  with  the  same  normal  trace  as  u{x)  on 

PROOF .     Let  u  be  the  solution  of 

'a{u,v)  =  f{v),       '^v  e'Hho.....hn 
i.e.,  a  is  given  by  equation  (3.11).  We  will  show  that 

From  this  follows  that 

a{u,v)  =  Y^a{u^,v}  =  a(u;  -\ \-  il,,v}  =  f{v),      Vv  £  Tin,. 

Since  this  holds  for  each  i  and  since 

"Who h,   =  'Hho  +  •  •  •  +  Hh,: 

follows  that 

-a{i,v)  =  f{v),       Went., ;,,. 

i.e.,  a{u,Tj)  is  represented  by  equation  (3.12). 

We  first  prove  equation  (3.13),  first,  for  /,  then  for  the  remaining  i.    Let  v  € 
y-ho hr  Then 

a{Piu,Piv)  =  a(P,tr,,P,t7)  +  a(P,(u,_i  +  •  •  •  +  tT;),  P/u). 
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Using  the  definition  of  {ui},  and  the  properties  of  orthogonal  projections,  we  obtain 

a{Piu,Piv)  =  a{{Li,v)  -f  a(u;_i  + \-ul,Piv). 

The  second  term  is  zero  since  u;_i  +  •  •  •  -f  u^  is  orthogonal  to  Tihr  Thus, 

a{Piu,Piv)  —  a{ui,v). 
Next,  we  consider  the  case  of  z  =  0, ...,/  —  1.  There  are  two  parts  to  prove. 

1.  Let 

Then,  since  the  support  of  v  is  in  Cli+i,  using  the  definition  gives  Eq^^v  =  0 
for  j  <  i.  Thus,  Mj{u,v)  =  0,  for  j  <  i. 

2.  Let 

We  split  u  in  four  parts, 

u  ^  {ui  + h  Ui+i)  4-  (tT,)  +  {u,_i  + \-u*)  +  (u*^i  + \-u*i). 

Since  the  first  and  fourth  parts  have  support  in  fii+i,  the  application  of  ^q^^, 
to  them  gives  zero.  Using  the  definition  of  {T/, . . . ,  tl,,  we  see  that 

Hence,  Hq      u,  =  u,  and  P,Hq^  ^u,  =  u,.  Similarly,  we  obtain  that 

a{ui_i  + \-u*,v)  =  0,      yveTih,, 

and 

^n.„  (^.-1  +  •••  +  «*)  =  {ui-i  +  •  •  •  +  ^). 
So  far,  we  have  obtained 

Mi{i,v)  =  a{u„PiHl^^^^v)  +  a{Pi{ui.i  +  •■■  +  u;),P,H},^^^v). 

But,  the  second  term  is  zero  since  PiH}^     v   G   "HKi,  and  u.-i  +  •  •  •  +  tT*  is 
orthogonal  to  Hh,-  Thus 
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Since  V  E  Hho  +  •  •  •  +  'Hh,,  it  follows  that 

which  is  orthogonal  to  iT,-.  Thus,  we  obtain 

Mi{u,v)  =  a{u,,{P,Hl^^^^v  -  v)  +  v)  =  a{ui,v).n 

We  now  present  the  definition  of  the  cosine  of  the  angle  between  spaces,  and  apply 
it  to  the  case  of  two  level  FAC  methods  to  get  a  bound  for  the  spectral  radius  of  the 
error  propagation  map.  See  Mandel  and  McCormick  [24]  and  Maitre  and  Musy  [23]. 

DEFINITION.     Let  U  and  V  denote  two  subspaces  of  a  Hilbert  space  7i  with  inner 
product  a(.,.).  We  define  the  cosine  of  the  angle  between  U  and  V  as: 

cos{u,v)=         sup 

u€U,vey.u.vj!iO  Ja{u,u)a[v,v) 

Note  that  cos{U,V)  G  [0,1],  and  if  U  n  V  is  nontrivial,  then  cos{U,V)  =  l.UU 
and  V  are  orthogonal,  then  cos{U.V)  =  0.  Furthermore,  if  tx  G  f/,  f  G  V,  then 


a(^l,^•)  <  cos{U,V)'^J a{u.  u)J a{v ,  v). 

When  cos{U,V)  <   1,  this  is  often  referred  to  as  a  strengthened  Cauchy  inequality. 
For  a  proof  of  the  following  two  Lemmas  see  Mandel  and  McCormick  [24]. 

Lemma  24  Let  U  and  V  be  nontnvial  subspaces  of  a  Hilbert  space  7i.   Then 

p{PuPv)  =  WPuPvPuW  =  cos^{U,V). 

where  Pu^Py  denote  the  orthogonal  projections  onto  the  respective  spaces,  and  p^PyPy) 
denotes  the  spectral  radius.    The  norm  is  the  Hilbert  space  operator  norm. 

Lemma  25   Let  U  and  V  be  nontnvial  subspaces  ofTl. 

1.  IfUnV  =  {0},  then  cos{U^,V^)  >  cos{U,V). 

2.  Ifn^  U  ®V,  then  cos{U^,V^)  =  cos{U,V). 

3.  IfU*  C  U  and  V*  C  V  and  H  =  U*  ®  V* ,  then  cos{U^,V^)  <  cos{U*,V*). 
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We  can  use  the  preceeding  two  Lemmas  to  estimate  the  spectral  radius  of  the  two 
level  FAC  method,  on  TihoM  using  the  subspaces  Tih^  and  7i/,i-  We  let 

U'  =  U  =  7^^^, 

and 

V*  =  {w  =  {i-  u'^)v  ■.venh,}cv  =  nh„ 

where  11''°  is  the  coarse  mesh  Raviart-Thomas  interpolation  map,  as  defined  in  Chap- 
ter I.  Thus,  the  spectral  radius 

To  obtain  a  bound  for  cos{7ikoi  [I  —  n''°)Whi))  ^^  work  on  each  coarse  mesh  element 
separately,  and  sum  over  all  the  coarse  mesh  elements.  Thus  suppose  that  for  each 
element  K  G  t^°{VIi),  there  exists  a  non-negative  constant 

0  <  <TA-   <   1, 

such  that 


Then 


aK{u,v)  <  cTK-^'aK[u.u\laK{v,v),    Vtx  €  Hh,,  and  Vr  e  (/  -  H'^")?^;.,, 


a{u,v)  =  ^aK{u,v)  <  Y^(TKyaK{u,u)^aK{v,v) 

K  K 


Thus, 


<  {rm.yiaK)\]a{u,u)yJa{v,v). 
cos{nHo,{I  -  n''°)7{,J  <  (max^K), 


K 

where  aj^  is  cos{Tiha\j^,{I  —  Il^°)7ih^\K)-  We  present  here  the  results  of  an  explicit 
computation  of  such  a  constant,  using  techniques  presented  in  Maitre  and  Musy  [23]. 
We  consider  the  particular  case  of  the  lowest  order  Raviart-Thomas  spaces,  i.e.,  r  —  0, 
defined  on  a  mesh  of  triangles. 

Let  K  denote  a  triangle  in  the  coarse  mesh  T^{Cli),  with  vertices  01,02,03.  Let  the 
angle  at  each  vertex  be  donoted  by  ^1,  ^2)^3)  respectively.  We  consider  the  refinement 
r'''(r2i)  obtained  by  dividing  each  coarse  triangle  K  into  four  equivalent  subtrian- 
gles  by  connecting  the  midpoint  of  each  edge  of  K.  We  denote  the  subtriangles  by 
Ki,. . .  ,K4.  See  Figure  3.2. 
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Figure  3.2:    Standard  refinement  of  a  triangle,  used  to  compute  the  angle  between 
coarse  and  fine  mesh  spaces. 


(0,0) 


Recall  that  a  function  in  Xh^iClo)  for  t-  =  0,  restricted  to  K,  is  specified  by  the 
three  values  of  the  normal  component  on  the  mid-point  of  edges.  But  for  functions  in 
Tiho,  ^h^  divergence  free  constraint  reduces  the  number  to  two.  Thus,  the  dimension 
of  TiholK  is  two.  A  basis  for  the  two  linearly  independent  divergence  free  functions  in 
TihoiK  is  given  by: 

ui  =  (1,0)  and  uj  =  (0,1). 

There  are  nine  edges  amongst  the  four  subtriangles,  and  there  are  four  constraints 
giving  us  five  linearly  independent  divergence  free  functions  defined  on  Ti^^  restricted 
to  K.  Since  two  of  them  belong  to  Ti^^,  there  are  only  three  basis  for  functions  in 
(/  -  Ii^°)nh,\K-  They  are  given  by: 

( 5 ,  — ^ 1       on  Ai 


U3 


U4 


■     /cos(g3)-cos(ga)      -Sm(g;)-5m(g3)\  ^^  J^ 

,  -cosfgjl-cosffl;)     -sin(g;)-l-sin(g3)\  ^^  J^ 

/cos(e])-cos(g;)     -sm(g;)-sin(tf3)  \  yy- 

/cos(^3)+cos(fl;)     5m[ejV-sin(«3j  \  „„  ^i' 
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and 

pcos(^g3)-l^sm[g^^      on   A'2 
(Cos(g3)-l^^imlg3l)      on   J^^ 

Computing  the  local  stiffness  matrix  on  if  using  the  basis  Ui^...^Ui^  we  obtain  a 
block  partitioned  matrix  of  the  form: 

Mk    Wk 
_Wl    Nk 

where 

'   (Ma-),j     =     aK[u,,Uj)\    i,j  =  1,2 

<    {Wk),j     =     aK[u,,Uj)]    t  =  l,2;    ;  =  3,4,5 
.   (-^A-)u      =     aA-(^.,'";);     2,i  =  3,4,5. 
On  A',  using  Lemma  24,  we  can  estimate  co5^(?Yho  U',  (^- H''")?^/,,  U)  by  the  spectral 
radius  of  the  error  propagation  map: 


co.'(Hh„|A-,(/-n''°)?^,jK) 


=  p\U- 


M^-'    0 
0  0 


Mk    Wk 

■T 
K 


n  Nk 


0    0 

0    Nji' 


Mk    Wk 
W^    Nk 


This  is  the  same  as 


p{M^'WkN^^W]:)  =  p{Ni:'Wl^MK'WK). 


In  the  case  of  the  Laplacian,  where 

A{x)  =  (    Q     ^    )    &nd  aK{u,v)  =  J    u'^vdx, 
we  obtain, 

s 


Mk  = 


4     0 
0    4 


■,Wk  =  j 


-1  -  C0S(^2)       C0S(^3)  -  C0S(^2)  COs(^3)  "  1 

-sin(^2)  -sin(^2) -sin(^3)     -  sin(^3) 


anc 


Nk  =  ^ 


2  1  1 
1  2  1 
1     1     2 


Here  S  denotes  the  area  of  triangle  K.  Furthermore, 

M^^'WkNk'W^ 
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1  +  cos2(^2)  +  cos2(^3)  sin(^2)  cos(^2)  -  sin(^3)  cos(^3) 

sin(^2)  cos(^2)  -  Ein(^3)  cos(^3)     sin'(^2)  +  sin^(^3) 


and  its  spectral  radius  is 


P  = 


+  ^|¥^Id[eM 


where 


1 
2 


d(^i,  ^2)  =  2  sin^(^2)  cos'(^3)  +  2  sin^(^3)  +  ;^  sin(2^2)  sin(2^3). 
For  an  equilateral  triangle 


and  for  a  right  triangle 


8' 
1 


P  = 


3.2.3      Quantitative  bounds  for  some  many  level  FAC  algo- 
rithms. 

Now  that  we  have  computed  the  spectral  radius  of  a  two  level  FAC  method,  we 
consider  an  application  to  the  case  of  the  many  level  FAC  method  using  the  same 
lowest  order  Raviart-Thomas  refined  mesh  spaces  based  on  triangles.  First,  we  present 
a  proof  of  Lemma  19,  for  the  FAC  algorithm.  Then,  using  this  Lemma  and  the  bound 
for  the  two  subspace  problems,  we  give  a  bound  independent  of  the  geometry. 

Recall  that  the  preconditioner  corresponding  to  the  symmetrised  FAC  algorithm 
is  given  by 

-a{u.,v)  =  aiPiv.Pa-^)  +  Y.a{PHl,^^^v,P,Hl,^^^v). 


t=i 


Note  that 


a(u,  v)  =  a{Piu,  Piu)  +  a{Pi^u,  Pi^u). 
We  can  form  a  preconditioner  by  replacing  P^^  by  H^^ .  We  obtain: 

-a,{u,v)  =  a{Piu,Piv)  +  a(i/^; V,  ^^7 V). 


Tl-l. 


Because  H^    u  is  the  ?^;_i  harmonic  extension  of  u  into  Cli,  and  because  Pi  u  —  Hq^  u 
is  a(.,  .)-orthogonal  to  Pi'^u.,  we  obtain  that 
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Thus, 

a2{u,u)  >  a{u,u). 

Note  that,  a2(., .)  is  the  preconditioner  corresponding  to  a  two  subspace  symmetrised 
multiplicative  Schwarz  algorithm  on  Hho,....h,  involving  the  subspaces  Tih,  and  X^'~o  Tih,- 
Thus,   we  would  require  solvers  to  determine  the  projections  onto  Ti^,   and  onto 

Ei=o  '^h,  ■ 

We  could  apply  the  same  technique  to  the  a{H'Q^^u,HQ^^u)  term,  i.e.,  decompose 

and  then  replace  Pi^-^Hq^^  by  H'q^I^Hq^^.    The  resulting  preconditioner,  which  we 
denote  by  asf., .)  is: 

a,{u,v)  =  a{Pu,Pv)  -f  a{P_,Hi,-'u,  Pi^^H'^-'v)  +  a{Hl,-'_^Hl,-'u,Hl,-'_Xy)- 

Note  that  H^^^l^Hl^^^  is  the  same  as  ^qT!,  •  Thus, 

az{u,u)  >  a2{u,u)  >  a{u,u). 

Similarly,  for  3  <  j  <  /  +  1,  we  define 

Note  that 

aui{u,v)  =  a{u,v). 

To  prove  Lemma  19,  we  consider: 

a{u,  u) 

_  ai+i{u,v)    ai{u,v)  d3{u,v)  a2{u^v) 

ai{u,u)    ai_i('u,ir)        a2{u,u)  a{u,u) 
Note  that,  by  construction, 

~     a{u,  u) 
To  find  an  upper  bound,  we  consider  each  term  in  the  product,  and  bound  it  by  using 

techniques  for  the  two  level  FAC  methods: 
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Since  only  the  last  terms  of  a_,+i(., .)  and  a_,(., .)  differ, 

Letting  w  =  Sq^'^    u  £  ?^o  +  •  •  •  +  "Wz-j+i,  we  need  a  bound  for 

a{Pi^,+^w,Pi-i+^w)  +  a{H^^^[^^w,H^^l^^w) 
a{w,'w) 

Thij   can  be  bounded  by  the  condition  number  of  the  two  level  symmetrised  FAC 

algorithm  applied  to  7Yo  +  •  •  •  +  '^^;-j+i,  with  the  subspaces  Tii-.j+i  and  T^oH VTi-i-j- 

This  completes  the  proof  of  Lemma  19,  in  the  case  of  iterative  refinement  methods. 
We  now  consider  a  quantitative  bound  for  the  many  level  symmetrisedYKC  algo- 
rithm for  the  case  of  the  lowest  order  Raviart-Thomas  finite  element  spaces  based  on 
triangles.  We  assume  that  the  refinements  in  the  appropriate  regions  are  obtained  by 
dividing  each  coarse  mesh  triangle  into  four  subtriangles  in  the  standard  way.  Then, 
each  of  the  two  subspace  symmetrised  FAC  has  a  condition  number  bounded  by 

1 

K.    <    . 

Thus,  the  condition  number  of  the  /  ^  1  level  symmttrised  FAC  is  bounded,  as  in 


Lemma  19,  bv 


1 


\-p\ 

For,  a  mesh  consisting  of  right  triangles,  the  spectral  radius  of  the  symmetrised  FAC 
method  is  /)  =  ^,  and  thus  k,  =  |,  thus  our  bound  for  the  /  -f  1  level  symmetrised 
FAC  method  in  this  case  is  given  by,  f|j   . 

For  the  case  of  a  mesh  containing  only  equilateral  triangles,  p  =  ^,  and  k  =  ||, 
and  the  condition  number  of  the  symmetrised  I  +  1  level  FAC  method  is  bounded  by 
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3.2.4      Step  3:   Determining  the  pressure. 

Let  {uh,ph)  denote  the  exact  solution  of  discrete  problem  (3.2).  Then,  by  construction, 

which  was  computed  in  the  previous  two  steps.  Recall  that,  for  the  mixed  formulation 
of  elliptic  Neumann  problems,  Wh  —  0,  and  so  ph  satisfies  the  following  equation: 

Find  ph  e  Yko,...,h,{^)  such  that 

<      a{ut,j  +  Uh,DF',v)  +  b{v,pH)     =    0  yveXh, h,{^)  (3.14) 

b{uh.i  +  Uh,DF,q)  =     Infi^M^)^^       Vg  e  iho, ...,/., • 

Now,  we  can  find  the  projection  onto  the  divergence  free  Raviart-Thomas  subspaces 

7ih,  by  substituting  the  new  residual  into  equation  (3.14).  For  i  =  0, . . . ,/  we  solve: 

Find  Wh.  e  XH,{n,),pH,  e  n.(a),  such  that 
a{wh,,v)  +  b{v,ph,)     =     -a{uh,i +  Uh,DF-,v)  =  b{v,ph),    VveA\(r2.)       (3.15) 

6(u\,g)  =  0,  vgGyH.(a). 

Each  u'h,  =  0,  but  Ph.  will  be  nonzero  in  general.  Using  the  fact  that  the  divergence 
operator  maps  A'';,,(fi,)  onto  F^,(n,),  we  see  that 


/. 


Ph  -  Phjqdx  =  0,      Vg  e  Yh,i^i), 
n, 

i.e.,  Ph,  is  the  i^  projection  of  the  discrete  pressure  solution  ph  onto  Yh,{Cli).  Recall 

that 

Yh,{ni)cL\n,)/Tom, 

implying  that  ph,  has  mean  value  zero  on  fi,-.  Thus,  if  we  knew  the  mean  value  oi  ph 
on  Qi,  then  we  would  know  the  L'  projection  of  p^  onto  Yh,{Cli)®  R,  where  R  denotes 
the  real  numbers.  For  each  i,  let  po^  denote  the  mean  value  of  ph  on  fi,-.  i.e., 

/n,  Phdx 

Jn.dx 

We  then  obtain 

'  {ph,  +  Pa  -  Ph)qd^  =  0,    Vg  e  n.(a)  e  'Po(a), 


i 


n, 

2 


i.e.,ph,  -\- Pq,  is  the  L    projection  of  p^  onto  Yh^{Qi)  ®  R.  Using  this,  we  can  compute 
the  mean  value  oi  ph  on  Q,+i: 

f  J         [      I        ,  \j  /n,+,(p/.. +Pn7M2; 

/       Phdx  =  /       [ph,  +pn,)dx  =*■  pn       = . 

•/n.+i  ^n.+,  Jn.^,  dx 
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Uiing  the  fact  that  the  discrete  pressures  are  discontinuous  across  element  boundaries, 
and  the  property  that  functions  in 

n„ H,{^)®Vo{n), 


w 


hen  restricted  to  r2,  —  i7,+i,  lie  in 


[n.(a)©7'o(a)]|n.-fi.,., 
we  obtain  that 

The  same  is  true  for  i   =   /,  when  f2/+i    =  0.  Thus,  if  ph,   and  pn7  are  known  for 
z  =  0, ...,/,  we  can  construct  ph' 

f   For  2  =  0,...,/ 

Ph  =  Ph,  +  Pn7,        on  Q,  -  fl,+i, 

-'n.  +  i 


and  pn.4, 


But.  for  i  =  0,  the  discrete  pressure  p^o  >  which  is  defined  on  the  whole  domain  Cl, 
has  the  same  mean  value  as  ph  on  Q,  namely,  zero.  Thus,  using  the  above  procedure, 
and  the  Z*  projections,  ph,,  we  can  construct  p/,. 


98 


Chapter  4 

A  Dirichlet-Neumann  algorithm 
for  Raviart-Thomas  elements. 


Let  fi  be  a  polygonal  domain  in  R'  which  is  triangulated  by  a  uniformly  shape  regular 
triangulation  t^{Q,).  Let 

f  n     =     fii  U  Qz  U  r,  where  Qi  fl  ^2  =  0, 

\  r   =5^1  n5r22, 

where  Vti  and  Q?  are  connected,  mutually  disjoint  subdomains.  We  assume  that  each 
fi,  is  a  union  of  elements  in  the  triangulation  T^{rt).  See  Figure  4.1. 

As  in  earlier  Chapters,  we  use  Xh{VL)  C  Ho,an{div,n)  and  Yh{Q.)  C  L^{Q)/Vo{^) 
to  denote  the  Raviart-Thomas  spaces  defined  on  the  triangulation  t''(Q),  for  an  eL 
liptic  problem  with  Neumann  boundary  conditions.  Recall  that  To{^)  denotes  the 
space  of  constants  on  Q.  We  decompose 

Xhin)  =  A'l  e  X2  e  A'3  ®  A'4, 

and 

Yh{n)  e  Vo{n)  =  Fi  ©  Fz  ©  1^3  ©  >4, 

where  the  subspaces  are  defined  by: 

Xi     =    Xh{n^)     =    XHin)nHo,anAdiv,Qi) 
X2    =    X,(n2)     =    Xh{n)nHo,anMv,^2) 

Y,    =  Ynin,)    =  YKin)nL'{£i,)/Vo{n,)  .^^^ 

Y2    =  YH{n2)    =  YH{n)nL^n2)/Vo{n2)  ^'' 

Y3    =  Vo{n,) 

Y,     =    ^(f^z). 


99 


Figure  4.1:  An  example  of  a  'T'  shaped  region. 
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Note  that  the  dimension  of  Y'3  as  well  as  of  ^4  is  one.  Before  we  define  X3  and  A''4  we 
define  Xh{T),  a  complementary  subspace  of  A'l  ®  X2  in  Xh{^)'- 

Xh{T)  =  {Hh  t  X^,{n)  :  Il\n,)u^,  =  0,  for  i  =  1,2}, 

where  n''(fi,)  denotes  the  interpolation  map  onto  Xh{^i)-    Let  tt^  be  a  function  in 
Xh{T)  which  has  a  normal  trace  of  non-zero  mean  value  on  F.  i.e., 


X 


ln'!^odS:r    r^   0- 


(4.2) 


For  instance,  we  could  define  tTq  to  have  a  normal  trace  of  1  on  F.  We  then  define 

X3  =  {cTfo  ■  c  E  R}. 
Note  that  the  dimension  of  A'3  is  1.  Next,  we  define 

A'4  =  {uh  e  A'h(F)  :   /  -TnUhds^  -  0}. 

Thus, 

AK(r)  =  A3®A4. 

Having  described  the  subspaces  A'l, .  . .  ,^4,  we  can  use  them  to  obtain  a  discreti- 
sation of  problem  (1.7).   By  choosing  a  basis  for  each  of  the  subspaces  and  ordering 
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them,  we  obtain  the  following  block  partitioned  stiffness  matrix: 


A„     0 
0        A22 

^13 


A^ 


A^ 


A'^ 

Bu  0 

0  B22 

0  0 

0  0 


Alz 
A23 
A33 

A^ 
^34 

-B13 

■B23 

JB33 

■B43 


A:4 
A24 
A34 
A44 

BlA 

B24 

0 

0 


BL    0 
0 

Biz 
Bfi    Bj^ 
0        0 


Biz 


0 
0 
0 


0 
0 
0 


0 
0 

^33 

0 
0 
0 
0 
0 


0 
0 

-^43 

0 
0 
0 
0 
0 


■  ui  ■ 

■  0   ■ 

'^2 

0 

Uz 

0 

Ui 

0 

Pi 

F, 

P2 

F2 

P3 

Fz 

.   P4    . 

F, 

(4.3) 


The  unknowns  in  Xi  are  represented  by  u,-,  and  those  in  Yi  are  represented  by  pi, 

in  the  matrix  notation.    We  note  that  A12  =  0,  since  functions  in  Xi  and  X2  have 

support  in  Qj  and  ^2  respectively,  which  are  mutually  disjoint.    Similarly,  Bi2,B2i 

are  zero  submatrices.  Next,  we  note  that  since  Yz  and  Y4  consist  of  piecewise  constant 

pressures,  and  since  the  velocities  in  A''i,A'^2)  and  X4  have  normal  traces  which  have 

mean  value  zero  on  each  dCl,,  we  obtain  that   Bzi,  B41,  Bz2i  B42,  Bzi,  B^  are  zero 

submatrices. 

REMARK.     Since  we  have  included  Vo{^)  (the  constants)  in  the  pressure  spaces,  the 

stiffness  matrix  will  have  a  null  space  of  dimension  one  corresponding  to  a  constant 

pressure. 

REMARK.     The  linear  system  (4.3),  has  the  form 


A     B^ 
B    0 


u 
P 


0 
F 


4.1      The  Glowinski- Wheeler  algorithm. 

In  this  section,  we  discuss  a  domain  decomposition  algorithm  due  to  Glowinski  and 
Wheeler  [16],  to  solve  equation  (4.3),  and  later,  in  the  next  section,  we  describe  a 
Dirichlet- Neumann  preconditioner.  We  describe  their  algorithm,  [16],  only  for  the 
case  of  two  non-overlapping  subdomains,  though  their  algorithm  is  applicable  to  the 
case  of  arbitrarily  many  non-overlapping  subdomains.  So  far,  the  Dirichlet-Neumann 
preconditioner  introduced  here,  holds  only  for  the  case  of  two  subdomains.  Each 
iteration  of  the  Glowinsk- Wheeler  algorithm,  involves  the  solution  of  a  saddle  point 
problem  on  each  of  the  two  subdomains. 

The  Glowinski-Wheeler  algorithm  is  divided  into  two  steps.    Let  (uh,Ph)  denote 
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the  solution  of  problem  (4.3).  We  decompose  the  discrete  solution  as 

where  /  denotes  initial  oi  inhomogeneous  and  H  denotes  pieciewise  discrete  harmonic 
and  divergence  free.  In  the  first  step  of  the  Glowinski-Wheeler  algorithm  we  deter- 
mine {uh,i,Ph.i)-  This  is  equivalent  to  the  saddle  point  problem  being  reduced  to 
a  symmetric  positive  definite  Schur  complement  problem  involving  the  unknowns  in 
X4.  i.e.,  the  velocity  unknowns  on  the  interface  F,  with  mean  value  zero  on  T.  In  the 
second  step  we  determine  {uh.H^Ph.n),  using  an  iterative  method,  such  as  the  con- 
jugate gradient  algorithm,  to  solve  the  positive  definite  Schur  complement  problem. 
Throughout  this  Chapter,  we  use  function  and  matrix  notation  interchangeably.  For 
instance,  we  use  the  vector  notation  [ui,  .  .  .,7)4]  to  denote  [uhiPh),  etc. 
To  describe  the  algorithm,  we  reorder  the  unknowns  as  follows: 

(rl) 
(r2) 
(r3) 
(r4) 
(r5) 
(r6) 
(r7) 
(r8) 


An 

BI, 

0 

0 

V4l3 

yli4 

0 

0 

Bn 

0 

0 

0 

5i3 

5:4 

0 

0 

0 

0 

A22 

Bf, 

A22 

A2A 

0 

0 

0 

0 

B22 

0 

B23 

^24 

0 

0 

5f3 

AL 

B22 

A33 

A34 

^3^3 

^4^3 

4^ 

Bf. 

Al, 

Bl. 

4^ 
-^34 

Aaa 

0 

0 

0 

0 

0 

0 

^33 

0 

0 

0 

0 

0 

0 

0 

543 

0 

0 

0 

'  ^1 

"  0   ■ 

Pi 

fl 

U2 

0 

P2 

F2 

^3 

0 

U4 

0 

P2 

F3 

.    P4    . 

F, 

(4.4) 


STEP  1  .  The  first  step  in  the  Glowinski-Wheeler  algorithm  is  to  find  an  initial 
discrete  velocity  and  pressure  {'Uhj,Ph.i)  G  A'h(n)  x  Yh{Vl),  such  that  when  this  is 
substituted  in  equation  (4.4),  and  the  residual  is  computed,  the  residual  is  zero  in  all 
rows  except  possibly  in  row  r6.  To  find  {uh,i,Ph.i)i  we  solve  a  series  of  problems  on 
substructures. 

To  start,  we  compute  Neumann  boundary  data  on  T  and  then  solve  local  problems 
on  ^1  and  ^2,  using  the  given  f{x)  in  the  divergence  constraint  and  the  computed 
Neumann  boundary  data  extended  by  zero  outside  T  as  the  specified  Neumann  bound- 
ary conditions  on  dQi  and  dCl2  respectively.  The  Neumann  boundary  data  must  be 
compatible  with  the  mean  value  of  f{x)  on  each  Q,.  After  the  two  local  problems  are 
solved,  we  compute  mean  values  for  the  pressures  in  Q.^  and  Q2  by  solving  a  system 
of  dimension  one,  which  thereby  reduces  the  residual  to  be  zero  with  respect  to  rb. 
We  outline  the  steps  in  more  detail: 
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step   la.        Define 

^.  ^    Jn,  f{x)dx 

Then,  c*7f^  will  provide  Neumann  boundary  data  on  F  which  is  compatible  with  f{x) 
on  both  fii,  since  by  assumption,  /(x)  has  mean  vcdue  zero  on  Q.  i.e., 

/    (/(x)  -  V  •  c'tTo)  <fx  =  0,      for  each  fi^. 

Step   lb.        Using  c'tTq,  we  solve  the  following  subproblems  for  z  =  1,2: 

Find  t7,,7  G  Xh(f],),  p,j  E  y/.(r2.)  such  that 

a{u,j,v)  +  b{v,pij)     =     -a(c*7ro,i;),  VveX^fi.)  (4.5) 

b{u,.i,q)  =     F(g)-fc(c-7fo,g),    Vg  G  ^(0.) 

By  the  construction  of  c'tTq,  the  local  problems  are  well  posed. 

Step   Ic.        We  choose  the  piecewise  constants  p^j  G  7^o(^i),  and  p^j  G  7^0(^2)  such 

that 

P3J  I     dx  +  p^j  I     dx  =  0,  (4.6) 

and 

a{c*£o^  uij +  U2,i,v)-r  b{v,p^j +P2J +  P3.I +  P4,i)  -  0,      {or  v  =  ifo-         (4.7) 
We  define 

UhJ   =   cVo  +  UlJ  +  U2,I, 

and 

P/i,/  =  Pl,7  +  P2,/  +  P3,/  +  P4,7- 

Note  that,  by  using  equation  (4.5),  and  the  fact  that  the  discrete  pressures  are  dis- 
continuous across  element  edges,  we  obtain  that 

BhUhj  =  Fh. 
We  now  give  a  block  Gaussian  elimination  interpretation  of  step  (1).  We  use 

[ui,I^PlJ.,U2,I,P2,I,U3j,U4j,p3j,Pij] 
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to  denote  the  matrix  representation  of  {uh,i,Ph,i),  obtained  after  the  initicd  step.  The 
equation  satisfied  by  {■Uhj,ph,i),  is: 


(rl) 
(r2) 
(r3) 
(r4) 
(r5) 
(r6) 
(r7) 
(r-8) 


An 
Bu 
0 
0 

^13 


Bf,  0 

0  0 

0  A22 

0  52: 
BL 


^23 

^r      dt      at 

"^14       -^^14       ^24 

0        0        0 
0        0        0 


0 
0 

-^22 


■^23 
-^24 


B\z    Bi4 
A23    A24 


B23 

A32 

A^ 
^34 

■B33 

■S43 


■B24 

>l34 
>l44 

0 
0 


0 
0 
0 
0 

Biz 
0 
0 
0 


0 
0 
0 
0 

-^43 

0 
0 
0 


Pi,/ 

■  0    ■ 

Fi 

U2,I 
P2,I 

0 
F2 

U3.I 

^4,7 

P3,7 

.  P4,7 

0 

F3 

F, 

(4.8) 


where  u^j  =  0. 

Finding  c*  in  step  (la)  is  equivalent  to  using  either  row  rl  or  rS  in  equation  (4.8) 
to  solve  for  u^j-  Equations  rl  and  rS  are  equivalent  since  the  matrix  has  a  null  space 
of  dimension  one  spanned  by  values  of  pa  and  p^  corresponding  to  a  constant  pressure 
on  Q.  Note  that  the  submatrices  ^33  and  B43  are  scalars,  and  nonzero.  Thus, 

B33    B43 

Solving  the  two  local  subproblems  in  step  (lb)  is  equivalent  to  solving  rows  rl,r2 
and  r3,r4  in  equation  (4.8)  for  Ui,/,?!,/  and  U2,/,p2,7  respectively,  by  substituting  in 
the  values  of  Uj  j  and  u^j  with  1/4  /  defined  to  be  zero,  i.e., 

-1 


p., 7 


A.    Bl 
B„    0 


-Ai3U3j 

Fi  -  B,2ll3.I 


(4.9) 


for  i  =  1,  2.  Thus,  after  steps  (la)  and  (lb),  the  equation  satisfied  by  {uh,l,Pi,l  +  P2j)i 
is: 

(rl) 

(r-3) 
(r4) 
(r5) 
(r6) 
(r7) 
(r8) 

because  u^j  is  defined  to  be  zero  and  psj  and  p^j  are  so  far  undetermined. 

Step   Ic.         Now,  we  note  that  scalars,  p^j  and  p^j  can  be  chosen  so  that  row  rb 

in  equation  (4.8)  is  satisfied  and  so  that  the  mean  value  of  the  pressure  ph,i  on  Cl  is 


Au 

BI^ 

0 

0 

A^3 

Al4 

0 

0 

Bu 

0 

0 

0 

B,3 

Bl4 

0 

0 

0 

0 

A22 

^2^2 

A23 

A24 

0 

0 

0 

0 

B22 

0 

B23 

B24 

0 

0 

^13 

Br3 
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■^23 

^2^3 

A33 

A34 

^3^3 

^4^3 

^14 

Bl 

A'^ 

^24 

Bl 

A^ 
^34 

A44 

0 

0 

0 

0 

0 

0 

■B33 

0 

0 

0 

0 

0 

0 

0 

B43 

0 

0 

0 

"  -"1,7  ' 

■  0 

Pl,7 

Fl 

U2.I 

0 

P2.I 

F2 

U3.I 

W3 

0 

w. 

0 

Fz 

0 

F4 
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zero.  This  involves  only  the  solution  of  two  equations,  because  the  block  matrices, 
B33  and  B43  are  non-zero  scalars,  and  the  requirement  that  phj  to  have  mean  value 
zero  is  equivalent  to  equation  (4.6),  i.e., 

I    ^3^3^3,7  +  -BJ3P4,/        =      -  ELML'^i.I  +  B^p.j]  -  A33U3J, 
\     \^l\P3.I  +  1^2  |P4,/       =       0. 

Since  the  last  two  columns  of  the  matrix  in  equation  (4.8)  is  nonzero  only  in  the 
fifth  row,  the  choice  of  p^j  and  p^j  will  not  alter  the  right  hand  side,  except  in  the 
fifth  row,  where,  by  the  choice  of  ^3,7  and  p4j  it  is  constructed  to  be  zero.  Thus  W4 
would  still  be  defined  by  the  use  of  row  r6  in  equation  (4.8),  i.e., 

W4  =  A'^^u^j  +  B^^pij  +  A^^U2,i  +  BJ4P2,/  +  Al^usj. 

Thus,  the  residual  after  step  (1),  is  given  by 

[0,0,0,0,0, -W,,0,Of.  (4.10) 

It  is  for  such  a  right  hand  side,  that  we  use  iterative  methods  in  the  second  step  of 
the  Glowinski-Wheeler  algorithm.  Define 

{uh,H,'Ph,H)  =  {^h  -  UhJ,Ph  -  Ph,l)- 

Then,  {uh,H,Ph,H)  satisfies: 


(-1) 
(r2) 

(r3) 

(r4) 

(r5) 

(r6) 

(r7) 

(r8) 


A^,  Bf,     0 

Bu  0 

0  0 

0  0 


A'^  B'^ 

■^13  -^IS 

4^  B^ 

^14  -°14 

0  0 

0  0 


0 
0        0 

-422       B22 

B22    0 

T        r>T 


A-'       B^ 

^23       •'-^23 


^24 

0 
0 


5J4 

0 
0 


^13       >ll4       0 

-B13    -B14    0 

^23       A2A       0 

■B23    B2A    0 

^33       ^34 
A44 

0 
0 


^34 
■B33 
-B43 


-S33 
0 

0 

0 


0 
0 
0 
0 

■^4^3 

0 
0 
0 


"  Ui^H 

■  0 

Pl,H 

0 

U2,H 

0 

P2,H 

0 

U3.H 

0 

P3.H 

0 

.  P4,H 

0 

(4.11) 


STEP  2 .         We  now  describe  a  Schur  complement  system  which  is  equivalent  to 

system  (4.11),  and  involves  the  unknowns  in  u^^jj. 

Step  2a.        By  using  rows  r7  or  r8  in  equation  (4.11),  we  determine  that 


U3,H 


0. 


Step  2b.        Following  that,  since  ^3^  =  0,  we  can  determine  Ui,H>Pi,fl'jt^2.H)P2,H  in 
terms  oiu^^u,  by  using  rows  rl,r2,r3,r4  in  equation  (4.11).  Namely, 


P>',H 


Bu    0 


-1  r 


Ai4 
Bi4 


U4H,      for  i  —  1,2. 
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Step  2c.  Since  u^fj  =  0,  once  iii,H,Pi.Hi'^2,HiP2.H  and  u^j}  are  known,  we  can 
determine  pa //  and  p^n  oi  mean  value  zero  on  Q  by  requiring  that  we  obtain  a  right 
hand  side  of  zero,  in  row  rb  of  equation  (4.11),  i.e.,  by  solving: 

\    |^l|P3,H  +   \^2\P4.H      =      0. 

Thus  far,  we  have  not  chosen  U4H  so  that  row  r6  of  equation  (4.11),  is  satisfied. 
Substituting  the  values  of  Ui^h,Pi,h-,'U2,HiP2.h  and  n^^fj  in  terms  of  ^4,^,  into  row 
r6  of  equation  (4.11),  we  see  that  u^^h  must  satisfy  the  following  Schur  complement 
system: 

Sru,,H  =  -W„  (4.12) 


where  5r  is  defined 

by 

5r  =  A44  — 

Ai4 
Bi4 

r 

_Bu 

Bl,] 
0 

-1 

■  >1:4  ■ 

B^4 

- 

A2i 
^24 

T 

A2_ 

B22    0 


Bj; 


A2i 

B24 

(4.13) 

Positive  definiteness  of   Sr-   We  now  show  that  ^r  is  positive  definite.   Let  U4 
be  arbitrarily  chosen,  and  let  ^''4  be  defined  by 

W4  =  SrU4. 

Let  Ui,pi,...,p4  be  constructed  using  U4  such  that: 


^1 

Bli 

0 

0 
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A^ 
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^3 
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A44 

0 

0 
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0 

0 

^4^3 

0 
0 

0 

0 

0 

0 

B43 

0 

0 

0 

'  ^\ 

■  0    ■ 

h 

0 

U2 

0 

P2 

0 

U3 

0 

U4 

W4 

h 

0 

.   P4    . 

0 

for  instance  by  using  steps  (2a),  (2b),  (2c).  Then  it  follows  that: 


-T  c  - 
•U4  orU4  = 


■  til  ■ 

r 

Pi 

U2 

P2 

^Z 

U4 

h 

.   P4    . 

^n 

B\, 

0 

0 

A,3 

Ai4 

0 

0 

5„ 

0 

0 

0 

B,3 

B\4 

0 

0 

0 

0 

A22 

^2^2 

A23 

A24 

0 

0 

0 

0 
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0 

B23 

B24 

0 

0 

A^ 

^13 

B13 

A^ 

"^23 

^2^3 

A33 

A34 

^3^3 

BL 

^14 

B14 

^24 

^2^4 

A'^ 
"^34 

A44 

0 
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0 

0 

0 

0 
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U4 

P3 
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(4.14) 


(4.15) 
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For  nonzero  U4,  the  discrete  velocity  [ui,U2,U3,U4\  in  equation  (4.15)  is  non-zero  and 
divergence  free.  Thus,  using  equation  (4.15)  and  the  fact  that 


if  B^u  =  0  then 


u 

P 


A    B^ 
B    0 


u 

V 


=  u  >lu, 


we  obtdn  that  u^SiUi^  >  0,  since  A  is  positve  definite  on  the  divergence  free  space 

H{div°,n). 

REMARK .     We  can  compute  the  action  of  Sr  on  any  vector  U4,  by  using  matrix  vector 

multiplications  and  by  solving  local  subproblems.    Thus,  we  can  use  the  conjugate 

gradient  algorithm  to  solve 

Sru,,H=-W„  (4.16) 

without  actually  computing  the  Schur  complement  Sr-  Note  that,  once  U4^h  is  known, 
the  rest  of  the  unknowns  can  be  computed  using  steps  (2a),  (2b),  and  (2c),  as  we  have 
already  described.  We  refer  to  this  basic  algorithm  as  the  Glowinski-Wheeler  algo- 
rithm. Results  of  numerical  tests  using  this  algorithm,  with  a  preconditioner,  shows 
that  the  rate  of  convergence  depends  on  the  mesh  width  /i,  cf.  Glowinski  and  Wheeler 
[16].  In  the  next  section,  we  describe  a  Dirichlei- Neumann  preconditioner,  which  we 
show,  leads  to  a  preconditioned  system  with  a  condition  number  independent  of  h. 

4.2      A  Dirichlet-Neumann  preconditioner. 

In  this  section  we  introduce  a  preconditioner  to  solve  system  (4.16)  in  the  second  step 
of  the  Glowinski-Wheeler  algorithm,  i.e.,  a  preconditioner  for  the  Schur  complement 
Sr  defined  by  equation  (4.13).  We  follow  ideas  described  in  Bj0rstad  and  Widlund  [4] 
and  Bramble,  Pasciak  and  Schatz  [5];  for  related  domain  decomposition  algorithms 
for  the  Stokes  problem,  see  Pasciak  [28j  and  Quarteroni  [29].  The  preconditioner  will 
involve  the  solution  of  a  saddle  point  problem  on  one  of  the  subdomains  Q;,  with 
mixed  Neumann  and  Dirichlet  boundary  conditions  on  the  boundary  5f2,-.  Therefore, 
we  first  discuss  the  saddle  point  formulation  of  elliptic  problems  with  mixed  boundary 
conditions. 

Consider  the  following  elliptic  problem  for  p  on  J^;: 

-V  •  {A{x)Vp)     =     0     on  a, 

n  ■  {A{x)Vp)    =    0    on  n  =  dn,  -  T, 
p    =    g   on  r,  c/  e  hII^{T). 
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The  saddle  point  formulation  of  this  problem  is: 

Find  u  G  ^o,r.^(<^i^',  ^.),  P  ^  L^i^i)  such  that 
a{u,v}  +  b{v,p)     =     J^gii-vds:,,    VtT  G  Ho,r=,  ((fii;,^,)  (4-17) 

b{u,q)  =  0,  \/qeL\n,), 

where  a{u,v)  =  /„  v^Aixy^udx,  and  6(v,p)  =  /^^pV  •  vdx,  as  in  the  previous  Chap- 
ters. Note  that  unlike  the  case  of  Neumann  boundary  conditions  on  the  whole  bound- 
ary dQ,i,  the  pressure  p  is  unique  in  Z^(Qi).  Also,  since  the  Neumann  data  is  prescribed 
only  on  T%  the  normal  trace  of  u  =  A{x)Vp  is  unknown  on  T.  This  will  be  reflected 
in  the  finite  element  discretisation  of  problem  (4.17). 

REMARK .  Recall  that  for  functions  in  H{div,  Q^)  the  normal  trace  is  defined  to  be  in 
H-'^'^{dQ.,),  by  Lemma  4.  Unfortunately,  the  restriction  of  a  H-'^l'^{dVti)  function  to 
a  subset  of  517,,  like  f  C  5Q,,  does  not  neccessarily  lie  in  /f"^/^(f ).  If  that  were  not 
the  case,  then  by  using  duality,  we  could  show  that  H'^I'^{T)  functions  extended  by 
zero  outside,  lie  in  H^^^{dQ,).  This,  we  know,  is  false. 

REMARK.  However,  it  can  be  shown  that  the  restriction  of  a  H~'^^^{dQ,)  function 
restricted  to  f  lies  in  (i7oo'(f))'.  This  is  because,  we  can  define  the  action  of  the 
restriction  on  any  function  in  Hqq  (F)  as  follows: 

where  <^  €  Hli\f),  and  E°(p  6  H^^^dQ,)  is  the  extension  by  zero  of  (f>.  Recall  that 
E°  is  a  bounded  map,  by  using  the  definition  of  /Too  (T).    See  Thomas  [32].   Thus, 
the  linear  functional  involving  g  in  equation  (4.17),  is  well  defined. 
REMARK.     For  any  f  C  5fi,,  the  space  H^fidiv^Q,)  is  defined  as: 

{ueH{div,n,)  :  {-f^u,4>)  =  0,    V<^e  ffoo'(f)}. 

The  discretisation  of  saddle  point  problem  (4.17)  is  obtained  by  restricting  the 

variational  problem  to  appropriate  Raviart-Thomas  finite  element  subspaces  oi Hox^^idiv,Qi) 

and  i2(Q.),  namely  T4(fi.)  n  Hoxiidiv^Q,)  and  Q;.(fi.)  n  L^Qi),  respectively.  Recall 

that  Fh(^i)  and  Qh{^i)  denote  the  Raviart-Thomas  spaces  for  the  Dirichlet  problem 

on  n,,  as  described  in  Chapter  1.  For  i  -  1,  we  obtain  the  foUowing  linear  system: 

Au  Bf,  A,z  Ai,  0 
Bn  0  5i3  5i4  0 
Aj,    B?3    Ag^       Ai\^    Bl         u,,H      =      0  (4.18) 


(r3) 
(r4) 
(r5) 


14       -^34 


^34 
^44 


0 


0 


B 


33 


■^33 

0 

0 


Ul,H 

■  0 

Pl,H 

0 

U3,H 

= 

0 

Ua.H 

G, 

P3.H 

0 
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For  i  =  2,  we  obtain: 


(rl): 

■    ^2:       -^22 

A23 

>l24 

0     • 

y-2.H 

■  0 

(r2): 

J522       0 

B23 

B24 

0 

P2,H 

0 

(r3): 

^23       -'-'23 

^33 

■^34 

Bl 

U3.H 

= 

0 

(r4): 

J^24      -^24 

.(2)7- 
^34 

/144 

0 

U4,B 

G4 

(r5): 

.  0        0 

■B33 

0 

0 

.  P^-^ 

0 

(4.19) 


We  have  used  the  notation  A^j    to  denote  the  appropriate  submatrix  in  the  stiffness 

matrix,  where  each  entry  is  obtained  by  integrating  the  appropriate  bilinear  form 

(2) 
only  over  subdomain  Qi.  Similarly  A\j'  is  obtained  by  integrating  over  ^2-  Thus, 

A,j  =  A,j   +  A^j  . 

Let  us  construct  a  Schur  complement  system  equivalent  to  system  (4.18),  by  using 
block  Gaussian  elimination  on  system  (4.18).  Using  row  rb  in  equation  (4.18),  we  find 
that  xi3H  —  0.  Next,  we  solve  for  the  rest  of  the  unknowns  in  terms  of  u^u-  Using 
rows  rl  and  r2  of  equation  (4.18),  we  obtain: 


Pl,H 


An     B^, 
Bn     0 


-1 


An 
Bii 


U4,H- 


Before  we  describe  the  equation  satisfied  by  u^e,  we  describe  how  the  mean  value  of 
the  pressure,  ps.Hi  is  determined  in  terms  of  the  rest  of  the  unknowns.  Substituting 
the  values  obtained  for  Ui^HiPi,Hi'^3M  ^^^  ^4,J7  i^^°  ^^'^  ^^  of  equation  (4.18),  we 
determine  psn  by 


P3.H  =  (-)- 


B^ 


33 


To  determine  the  equation  satisfied  by  1x4,^5  we  substitute  the  values  oiui^jj,pi^ij 
and  u^H,  in  terms  of  u^jj,  into  row  r4  of  equation  (4.18).  We  then  obtain  that 


where 


S'r 


(1) 


A 


(1) 
44 


Al4 
B\4 


■[T    r 


^11     0 


1  -1  r 


AiA 
Bi4 


:-(i) 


(1) 


Thus  5p     is  the  Schur  complement  associated  with  u^^h  in  system  (4.18).    5p     is 


.rc(i 


positive  definite  since  u^S^  u^  equals  the  a(.,.)  norm  of  the  divergence  free  discrete 
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velocity  [ui.h, '"s.i/,  ^4,h]  obtained  by  solving  system  (4.18)  and  that  is  positive  defi- 
nite; the  proof  is  identical  to  that  used  to  show  that  5r  is  positive  definite. 
Similarly,  using  block  Gaussian  elimination  on  system  (4.19),  we  obtain 


Sr  '^A.H  =  G4, 


where 


:(2) 


A 


(2) 
44 


A24 
B24 


A22       B22 

B22    0 


A2A 

-B24 


is  the  Schur  complement  associated  with  Uj^jj  for  the  mixed  boundary  condition  prob- 
lem  on  subdomain  Qt,  namely  problem  (4.19).  5p     is  also  positive  definite. 
Since, 


A..  -  A^^^  +  A^^^ 

j^44   —   ^44     "I-  .(^44  , 


it  follows  that 


Sr  =  S- 


(1) 


s''\ 


(4.20) 


:(^) 


Thus,  we  could  use  either  S^'  or  5p  as  a  preconditioner  for  the  Schur  complement 
Sy-  Note  that,  solving  linear  systems  involving  5p  or  5^  ,  can  be  done  without 
computing  Sj-  or  S^  ,  by  solving  systems  of  the  form  (4.18)  or  (4.19),  respectively. 
We  now  consider  the  condition  number  of  the  preconditioned  system.  This  deter- 
mines the  rate  of  convergence  of  the  preconditioned  conjugate  gradient  algorithm. 

Theorem  7  Fori  =  1,2,  consider  the  preconditioner  S-p    for  the  matrix  Sr-    There 
exists  a  positive  constant  Ci{Qi,Cl2)  independent  of  h  such  that 


1  < 


u'^fjSrii^^H 


<  1  +  c,. 


(4.21) 


PROOF.      To  be  specific,  let  i   =    1.     A  similar  argument  holds  for  i   =   2.     Using 
equation  (4.20),  we  see  that 


T      c(l) 


(1), 


>  1. 


u\jjS\^  'Ui^H 


This  proves  the  lower  bound. 

To  obtain  the  upper  bound,  we  consider 


(2), 


^4,H'-'r  ^4,H 


=  1  + 


(1) 


T       P\^ ) 
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Recall  that  u^ ^Stu^^h  is  the  a(., .)  norm  of  the  divergence  free  discrete  velocity  exten- 
sion [ui^Hi U2.niU3,Hi U4,h]i  which  is  the  velocity  part  of  the  solution  of  equation  (4.11) 
with  —W4  defined  by  —W4  =  SrU4_H-  Similarly,  u'^jjSl  u^^jj  is  the  a(.,.)  norm  of  the 
divergence  free  discrete  velocity  extension  [tii,if)^i3,F?^'4,H])  which  is  the  velocity  part 
of  the  solution  of  equation  (4.18)  with  G4  defined  by  S^  u^^h.  u'^gSj'  u^^h  is  the 
a(.,.)  norm  of  the  divergence  free  discrete  velocity  extension  [u2,Hi^3,H,^i4,H],  which 
is  the  velocity  part  of  the  solution  of  equation  (4.19)  with  G4  defined  by  5^  'u/^^h- 
We  seek  to  bound 

for  some  positive  constant  Ci{Q,i,^2).  Since  row  rl  of  equation  (4.19)  is  homogeneous, 
it  follows  that  [u2,/f,'"3,H5  ^^4,i/]  is  a(.,  .)-orthogonal  to  all  divergence  free  discrete  ve- 
locities  in  Xh{^2)  =  1^h(f^2)  l^-^o.an2(^"')  ^2)-  Thus,  u\hSv  U4,h  is  less  than  or  equal 
to  the  a(.,.)  norm  of  any  discrete  divergence  free  function  in  Vh{Q2)  with  the  same 
normal  trace  as  [u2,H5  ^3,H)  ^4,//]-  Note  that,  the  normal  trace  of  [ti2,H»'^3,lf  )^4,if] 
is  zero  on  Fj  =  9^2  —  T,  and  equals  u^h  on  T  (up  to  sign).  Thus,  in  particular, 
uj jjSy  i^4,H  is  less  than  or  equal  to  the  a(.,.)  norm  of  the  discrete  divergence  free 
extension  f?''(7„[u2,H,"3,H, ti4,H])|an25  given  by  theorem  2  (the  Extension  theorem), 
in  Chapter  1.  By  the  Extension  theorem,  the  a(., .)  norm  of  the  extension  is  bounded 
by  the  H''^^^{dQ,2)  norm  of  the  normal  trace,  with  a  constant  independent  of  h. 

Next,  we  note  that  the  H~'^^'^{dQ2)  norm  of  7„[u2,if,^i3,//,  ■"4,H]|sn2  can  be  bounded 
by  the  H'''^^^{dQi)  norm  of  7n[^ii,H,'"2,H,^i3,H5«4,H]lsni5  with  a  constant  independent 
of  /i,  but  possibly  depending  on  the  two  subdomains.  This  is  because,  we  can  map 
Sr^i  onto  00,2  by  a  bijective  mapping,  keeping  T  fixed.  This  establishes  an  equivalence 
between  H'^^^{dQ,i)  and  H^/'^{dO,2).  By  using  duality,  we  obtain  the  result.  Next,  we 
note  that  we  can  bound  the  H~^^^{dQ,i)  norm  of  7„[ui,H>'"2,H>^3,Hi«4,H]|ani  by  the 
a(., .)  norm  of  [ui^h,U2^h^'^4,h]i  with  a  constant  independent  of  h,  by  using  Lemma  4 
(the  normal  trace  Lemma).  But  the  square  of  this  norm  is  equal  to  u^jjSj'  U4,H- 
Thus,  there  exists  a  positive  constant  c^(Qi,f22),  independent  of  h,  such  that  the  up- 
per bound  in  equation  (4.21)  holds.  □ 

REMARK.  Note  that  Sy\  maps  the  Neumann  data  on  T  onto  Dirichlet  data  on  T. 
This  is  done  by  finding  the  Dirichlet  data  of  a  discrete  velocity  satisfying  the  piece- 
wise  homogeneous  equation  in  (4.18)  and  satisfying  the  given  Neumann  boundary 
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condition  on  T.  Thus,  we  may  refer  to  5p    as  a  Neumann-Dirichlet  map.  Since 

we  may  expect  heuristically  that  its  condition  number  be  O(^),  since  the  difference 
in  exponents  of  the  spaces  is  one,  corresponding  to  a  derivitive.  Since  Sr  is  spec- 
trally equivalent  to  5}-  ,  we  also  expect  Sr  to  have  a  condition  number  of  O(^).  Since 
5r  =  5r  +5p  ,  we  see  that  Sr  maps  the  Neumann  data  on  T  to  the  jump  in  Dirichlet 
data  on  F,  obtained  by  using  the  two  Neumann-Dirichlet  maps  on  the  two  subdo- 
mains.  Thus,  system  (4.12)  has  as  its  solution,  the  appropriate  Neumann  data  that 
leads  to  the  correct  jump  in  the  Dirichlet  data.  For  details  about  the  terminology  of 
the  Dirichlet-Neumann  algorithm,  see  Bj0rstad  and  Widlund  [4]. 
REMARK  ON  QUANTITATIVE  BOUNDS  .  We  note  that  we  can  easily  obtain  some  quan- 
titative bounds  for  the  condition  number  of  the  Dirichlet- Neumann  algorithm  in  the 
case  of  special  geometries  and  special  coefficients  for  the  elliptic  operator.  Suppose 
that  ^2  is  the  mirror  image  of  Qi,  under  reflection  across  F.  Further,  assume  that 
the  mesh  on  Q2  aJso  is  a  reflection  of  the  mesh  on  fli,  and  that  the  coefficients  A{x) 
are  even  across  F.  Then,  we  can  bound  u'^ fj^r  "4,if  by  u^ h^v  "4,h  since  the  even 
extension  of  [ui^h,'U-3,h^'^a,h]  across  F  into  ^2  provides  a  divergence  free  harmonic 
extension  as  used  in  the  proof  of  Theorem  7.  This  shows  the  condition  number  to  be 
bounded  by  2.  The  same  bound  holds  if  the  mirror  image  across  F,  of  the  mesh  on 
fij  is  contained  in  the  mesh  on  Qt- 
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